Jamming II: A phase diagram for jammed matter 



00 

o 
o 

(N 
< 

in 



O 



o 
o 



> 

T — I 
(N 

00 

O 

00 

o 



X 



Ping Wang^, Chaoming Song^ and Hernan A. Makse^'^ 
^ Levich Institute and Ptiysics Department, City College of New York, New Yorit, NY 10031, US 
^ Departamento de Fisica, Universidade Federal do Ceard, 60451-970 Fortaleza, Ceard, Brazil 

The problem of finding tlie most efficient way to pacic splieres has an illustrious 
history, dating back to the crystalline arrays conjectured by Kepler and the random 
geometries explored by Bernal in the 60's. This problem finds applications spanning 
from the mathematician's pencil, the processing of granular materials, the jamming 
and glass transitions, all the way to fruit packing in every grocery. There are presently 
numerous experiments showing that the loosest way to pack spheres gives a density of 
~ 55% (named random loose packing, RLP) while filling all the loose voids results in 
a maximum density of ~ 63 — 64% (named random close packing, RCP). While those 
values seem robustly true, to this date there is no physical explanation or theoreti- 
cal prediction for them. Here we show that random packings of monodisperse hard 



spheres in 3d can pack between the densities 



or 53.6% and 



or 63.4%, defin- 
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ing RLP and RCP, respectively. The reason for these limits arises from a statistical 
picture of jammed states in which the RCP can be interpreted as the ground state of 
the ensemble of jammed matter with zero compactivity, while the RLP arises in the 
infinite compactivity limit. We combine an extended statistical mechanics approach 
'a la Edwards' (where the role traditionally played by the energy and temperature in 
thermal systems is substituted by the volume and compactivity) with a constraint on 
mechanical stability imposed by the isostatic condition. We show how such approaches 
can bring results that can be compared to experiments and allow for an exploitation 
of the statistical mechanics framework. The key result is the use of a relation between 
the local Voronoi volumes of the constituent grains (denoted the volume function) and 
the number of neighbors in contact that permits to simple combine the two approaches 
to develop a theory of volume fluctuations in jammed matter. Ultimately, our results 
lead to a phase diagram that provides a unifying view of the disordered hard sphere 
packing problem and further shedding light on a diverse spectrum of data, including 
the RLP state. Theoretical results are well reproduced by numerical simulations that 
confirm the essential role played by friction in determining both the RLP and RCP 
limits. The RLP values depend on friction, explaining why varied experimental results 
can be obtained. 



I. INTRODUCTION 

For many years, a great deal of research has been de- 
voted to the properties of packings of granular materials 
and other jammed systems such as compressed emulsions, 
dense colloids and glasses jll, [1]. Generally speaking, 
there are three interrelated approaches to understand the 
nature of the jammed state of matter and its ensuing 
jamming transition: 

(a) Statistical mechanics of jammed matter: 
The jammed phase is described with the ensemble of 
volume and force fluctuations proposed by Edwards 
@, i, i, i, 0, i, i, 0, El • Ideas from the physics of glasses 
provide interesting cross-fertilization between the jam- 
ming transition and the glass transition [l^, [H, ■ 

(b) Critical phenomena of deformable particle: 
jamming is seen as a critical point at which observables 
such as pressure, coordination number and elastic mod- 
uli behave as power-law near the critical volume fraction, 
^c, Slllllll,!!! [11]. The jammed state 
is modelled as granular media composed of soft particles, 
typically Hertz-Mindlin ,22, 25, 26] grains under external 
pressure and emulsions compressed under osmotic pres- 
sure [131, as they approach the jamming transition from 



the solid phase above the critical density: 0^0+. 

(c) Hard-sphere glasses: In parallel to studies in the 
field of jamming there exists a community attempting to 
understand the packi ng pro blem approaching jamming 
from the liquid phase (28l. l29l [30l| . that is, — + Here, 
amorphous jammed packings are seen as infinite pressure 
glassy states [H, |31|. Therefore, the properties of the 
jamming transition are intimately related to those of the 
glass transition [2^. 

A great deal of effort has been devoted to the study 
of jamming applying these approaches. A large body of 
experiments and simulations have fully characterized the 
jammed state in the critical phenomena framework of 
soft-particles. The hard-sphere field studies frictionless 
amorphous packings interacting with hard-core normal 
forces only. Their results are a particular case of the 
more general problem of jammed matter composed of 
frictional and frictionless hard spheres as treated in the 
present work. 

Statistical mechanics studies of jammed matter were 
initiated in 1989 with the work of Sir Sam Edwards. It 
was first recognized that the main theoretical difficulty 
to develop a statistical formulation results from the lack 
of well-defined conservation laws on which an ensemble 
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description of jammed matter could be based. Unlike 
equilibrium statistical mechanics where energy conser- 
vation holds, granular matter dissipates energy through 
frictional inter-particle forces. Therefore, it is doubtful 
that energy in granular matter could describe the mi- 
crostates of the system and a new ensemble needs to be 
considered. 

Stemming from the fact that it is possible to explore 
different jammed configuration at a given volume in sys- 
tematic experiments [a H, [13, [13] , Edwards proposed the 
volume ensemble (V-ensemble) as a replacement of the 
energy ensemble in equilibrium systems. A simple exper- 
iment is merely pouring grains in a fixed volume and ap- 
plying perturbations, sound waves or gentle tapping, to 
explore the configurational changes. Concomitant with a 
given network geometry there is a distribution of contact 
forces or stresses in the particulate medium. This means 
that the V-ensemble must be supplemented by the force 
or stress ensemble (F-ensemble) determined by the con- 
tact forces [l,[3^[3j for a full characterization of jammed 
matter. 

Following this theoretical framework, a large body of 
experiments, theory and simulations have focused on the 
study of volume fiuctuations in granular media. Statis- 
tical studies have been concerned with testing for the 
existence of thermodynamics quantities such as effective 
temperatures and compactivity as well as challenging the 
foundations based on the ergodic hypothesis or equal 
probability of the jammed states. 

Unfortunately, there is no first principle derivation of 
the granular statistical mechanics analogous to the Li- 
ouville theorem for equilibrium systems [355 . Therefore, 
advancing the statistical mechanics of granular matter 
requires well-defined theoretical predictions that can be 
tested experimentally or numerically. While the possi- 
bility of a thermodynamic principle describing jammed 
matter is recognized as a sensible line of research, the 
problem with the statistical approach is that after al- 
most 20 years there are no practical applications yet. We 
require predictions of practical importance such as equa- 
tions of state relating the observables: pressure, volume, 
coordination number, entropy, etc, that may lead to new 
phenomena to be discovered and tested experimentally 
that will allow for a concrete exploitation of the thermo- 
dynamic framework. 

Here, we explore this problem by developing a sys- 
tematic study of the V-ensemble of jammed matter. Our 
systems of interest are primarily packings of granular ma- 
terials, frictional colloids, infinitely rough grains and fric- 
tionless droplets mimicking concentrated emulsions. We 
attempt to answer the following questions using statisti- 
cal mechanics methods: 

— What is a jammed state? 

— What are the proper state variables to describe the 
granular system at the ensemble level? 

— What are the ensuing equations of state to relate 
the different observables? For instance, the simplest 
equilibrium thermal system, the ideal gas, is described 



by pressure, volume and temperature {p, V, T) through 
pV = nRT. Is it possible to identify the state variables 
for jammed matter which can describe the system in a 
specified region of the phase space? 

— Which are the states with the maximum and mini- 
mum entropy (most and least disordered)? 

— Can we define a " compactimeter" to measure com- 
pactivity? 

— Can we use the ensemble formalism to define the 
random close packing and random loose packing state 
and predict their density value? 

Of similar importance is the study of how uniform 
spheres will pack if assembled randomly. This problem 
is the main application of the theory presented in the 
present study, and is discussed next. 



A. Sphere packing problem 

Filling containers with balls is one of the oldest math- 
ematical puzzles known to scientists. The study of this 
problem started four centuries ago when Johannes Ke- 
pler conjectured that the most efficient arrangement of 
spheres is the FCC lattice (an important part of the 18th 
problem proposed by Hilbert in 1900). Even though this 
is a tool used for centuries in fruit markets around the 
globe, nearly 400 years passed before this conjecture was 
considered a mathematical proof, which has been devel- 
oped only recently by Hales [3a| in a series of articles 
covering 250 pages supplemented by 3Gb of computer 
code to determine the best ordered packing through lin- 
ear programming. The difficulty arises since in 3d it is 
not enough to look at the packing of one cell, but it is 
necessary to consider several Voronoi cells at once. That 
is, the packing that minimize the volume locally (the do- 
decahedron) does not tile the system globally. Such a 
situation does not arise in 2d, where the hexagonal pack- 
ing minimizes the volume locally and globally; the Kepler 
conjecture in 2d was proved long ago. 

The analogous problem for disordered packings has 
also an illustrious but unfinished history. This prob- 
lem was initiated by the pioneering work of Bernal in 
the 1960's [33| (although earher attempts can be found). 
The traditional view states that [11] "packings of spher- 
ical particles have been shaken, settled in different flu- 
ids and kneaded inside rubber balloons and all with 
no better results than a maximum density of 63%" . 
This is the so-called random close packing RCP limit 
[13, [H, [H, [13, [H. On the other hand, other experi- 
ments have shown that densities as low as 55% can be 
obtained in random loose packings, RLP [sol [4ll. [i^ . 

To this date there is no physical explanation of this 
phenomenon, no theoretical prediction of such density 
values and heated debates are still found in the literature 
regarding the existence of rigorous definitions of the RCP 
and RLP, the uniqueness of the RCP state and the nature 
of their state of randomness. 
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B. Objectives 

Here we apply Edwards statistical mechanics to study 
the situation in random disordered systems. We argue 
that part of the solution will be found by treating these 
systems at the ensemble level. We illuminate a diverse 
spectrum of data on sphere packings through a statisti- 
cal theory of disordered jammed systems, while revealing 
new packing states open to experimentation. Our results 
ultimately lead to a phase diagram providing a unify- 
ing view of the disordered sphere packing problem that: 

( i ) provides a statistical interpretation of RLP and RCP, 

(ii) predicts their density values in good agreement with 
experiments, and (in) establish the concomitant equa- 
tions of state relating observables such as the average 
mechanical coordination number Z , entropy S*, and vol- 
ume fraction 0. 

The RCP limit arises as the ground state (zero com- 
pactivity) of jammed hard spheres predicted to be at 



the understanding of how random packings fill space in 
3d for any friction, and eventually can be extended as a 
mathematical model for packings in 2d and in higher di- 
mensions, anisotropic particles and other systems. While 
physical application may be limited to 2d or 3d, higher 
dimensions find application in the science of information 
theory, where data transmission and encryption find use 
for packing optimization. 

The packing problem is by its very nature an interdis- 
ciplinary subject with many important applications in a 
variety of systems, including liquids, amorphous solids, 
granular media, colloids, glasses and living cells. There- 
fore, the topic treated here finds interest in its own right 
in diverse areas from physics, chemistry, mathematics to 
biology. 



C. Outline 
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0.634, 



(1) 



with different Z. Thus, the RCP state is not a unique 
point in the phase diagram (0, Z), but there exist many 
RCP with different coordination number determined by 
the interparticle friction coefficient but with the same 
volume fraction. 

The RLP states appear in the phase diagram along a 
line of infinite compactivity parameterized by Z, ending 
at the minimum density predicted to be 
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0.536. 



(2) 



The nature of the disorder of the packings is statisti- 
cally characterized by the entropy, which is shown to be 
larger in the random loose case than in the random close 
case. 

The fundamental results that allow to describe the 
RCP and RLP states have been obtained in the compan- 
ion paper Jamming I [43| and describes jammed matter 
through a volume function, amenable to analytical cal- 
culations. Using this volume function, we combine an 
extended statistical mechanics approach 'a la Edwards' 
with a constraint on mechanical stability imposed by the 
isostatic condition. The key idea is the use of a rela- 
tion between the local Voronoi volumes and the number 
of neighbors in contact that permits to simply combine 
the two approaches. The analytical expression of the free 
volume function shows the use of the compactivity to de- 
scribe granular packings, a rather abstract concept up to 
now. 

An extensive program of computer simulations of fric- 
tional and frictionless granular media tests the predic- 
tions and assumptions of the theory, finding good general 
agreement. The phase diagram introduced here allows 



This paper is the second installment of a series of pa- 
pers dealing with different aspects of jammed matter 
from a statistics mechanics point of view. In a compan- 
ion paper. Jamming I [i^, we describe the definition of 
the volume function based on the Voronoi volume of a 
particle and its relation to the coordination number. 

In the present paper, Jamming II, we describe calcu- 
lations leading to the phase diagram of jammed matter 
and numerical simulations to test the theoretical results; 
a short version of the present paper has been recently 
published in [31 . Jamming III [455 describes the en- 
tropy calculation to characterize randomness in jammed 
matter and a discussion of the V-F ensemble. Jam- 
ming IV [i^ studies the distribution of volumes from 
the mesoscopic ensemble point of view and a generalized 
z-ensemble to understand the distribution of coordina- 



tion number. Jamming V 47[ details the calculations for 
jamming in two dimensions at the mesoscopic level. Fi- 
nally, Jamming VI [48] attempts a solution of the volume 
distribution by brute force, exact up to the second coor- 
dination shell of particles in disordered packings, which 
may overcome many difficulties that we will encounter in 
the present mesoscopic approach. 

The outline of the present paper is as follows: Section 
im explain the basis of the Edwards statistical mechanics. 
Section [ml summarizes the results of Jamming I regard- 
ing the calculation of the volume function. Section HVl 
describes the isostatic condition that defines the ensem- 
ble of jammed matter through the 6jam function, and 
explains the difference between the geometrical, z, and 
mechanical, Z , coordination number, which is important 
to define the canonical partition function. Section |V] de- 
fines the density of states in the partition function. Sec- 
tion IVII explains the main results of this paper in the 
equations of state and and the phase diagram 
of Fig. 21 Finally, Section IVIII explains the numerical 
studies to test the theoretical predictions. 
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II. STATISTICAL MECHANICS OF JAMMED 
MATTER 

Conventional statistical mechanics uses the ergodic hy- 
pothesis to derive the microcanonical and canonical en- 
sembles, based on the quantities conserved, typically the 
energy E [35| . Thus the entropy in the microcanonical 
ensemble is S{E) = ks log / 5{E — H{p, q))dpdq, where 
Ti-ip, q) is the Hamiltonian. This becomes the canonical 
ensemble with exp[—H{dS/dE)]. 

The analogous development of a statistical mechanics 
of granular and other jammed materials presents many 
difhculties. Firstly, the macroscopic size of the consti- 
tutive particles forbids the equilibrium thermalization of 
the system. Secondly, the fact that energy is constantly 
dissipated via frictional interparticle forces further ren- 
ders the problem outside the realm of equilibrium statis- 
tical mechanics due to the lack of energy conservation. In 
the absence of energy conservation laws, a new ensemble 
is needed in order to describe the system properties. 

Following this theoretical perspective, Edwards pro- 
posed the statistical mechanics of jammed matter and 
interest in the problem of volume fluctuations has flour- 
ished Q . The central concept is that of a volume function 
W replacing the role of the Hamiltonian in describing 
the microstates of the system in the V-ensemble and the 
stress boundary 



n, 



aij dV 



with the stress an 
ensemble (gI IssI. a4 



= 1/(21^) Ec/f^l describing the F- 
where f^yrf are the force and posi- 
tion at contact c. For simplicity only the isotropic case is 
described. Thus, only the pressure a = 0-^^/3 is necessary 
to describe force fluctuations. 

If we partition the space associating each particle with 
its surrounding volume, Wj (for instance with a Voronoi 
tessellation as it will be done bellow), then the total vol- 
ume, W, of a system of N particles is given by: 



N 



(3) 



The ensemble average of the volume function W provides 
the volume of the system, V = (W) , in an analogous way 
to the average of the Hamiltonian is the energy in the 
canonical ensemble of equilibrium statistical mechanics. 

The full canonical partition function in the V-F ensem- 
ble is the starting point of the statistical analysis 0: 



Z{X,A) ^ / .g(>V,H) exp 



n 

A 



X 



Gjam dW dn, 



(4) 

where g(yV, H) is the density of states for a given volume 
and boundary stress. Here Gjam formally imposes the 
jamming restriction and therefore defines the ensemble of 
jammed matter. This crucial function will be discussed 



at length below. As a minimum requirement it should 
ensure touching grains, and obedience to Newton's force 
laws. 

Just as dE/dS ^ T is the temperature in equilibrium 
systems, the temperature-like variables in granular sys- 
tems are the compactivity 



X 



dV 

as' 



(5) 



and the angoricity [from the Greek "ajx^C' (ankhos) 
stress! M: 



_ dll 



(6) 



[Note that in the angoricity is denoted Z. Here we 
use A for the isotropic angoricity since Z is used to de- 
note the mechanical coordination number. In general the 
angoricity is a tensor, Aij = here we simplify to the 
isotropic case]. 

The compactivity measures the power to compactify 
while the angoricity measures the power to stress. These 
quantities have remained quite abstract so far, perhaps, 
this fact being the primary reason for the great deal 
of controversy surrounding the statistical mechanics of 
jammed matter. In the present paper we attempt to pro- 
vide a meaning and interpretation for the compactivity 
(the angoricity will be treated in subsequent papers) by 
developing equations of states as well interpretations in 
terms of a "compactimeter" to measure the "tempera- 
ture" of jammed matter. In Eq. fl]), the analogue of 
the Boltzmann constant are set to unity for simplicity, 
implying that the compactivity is measured in units of 
volume and the angoricity in units of boundary stress 
(stress times volume). 

In the limit of vanishing angoricity, A 0, the sys- 
tem is described by the V-ensemble alone. This is the 
hard sphere limit which will be the focus of the present 
work where the following partition function describes the 
statistics: 



dW, 



(7) 



where g{W) is the density of states for a given volume 
W. 

From Eq. ([7]) we identify three minimal steps in devel- 
oping analytical solutions as discussed below. Section lllll 
discusses the need for a volume function, W. Section HVl 
discusses the need for a proper definition of the ensemble 
of jammed state allowing a definition of the constraint 
function Ojam- Section |V] discusses the derivation of the 
density of states. 

It is recognized that thermodynamic analogies may il- 
luminate methods for attempting to solve certain gran- 
ular problems, but will inevitably fail at some point in 
their application. The mode of this failure is an interest- 
ing phenomenon, as the range of phenomena for which 
ergodicity of some kind or other will apply or not is an 



interesting question. This has been illustrated by the 
compaction experiments of the groups of Chicago, Texas, 
Rennes and Schlumberger [1, [1, [13, [SI] . They have shown 
that reversible states exist along a branch of compaction 
curve where statistical mechanics is more likely to work. 
Conversely, experiments also showed a branch of irre- 
versibility where the statistical framework is not expected 
to work. Poorly consolidated formations, such as a sand- 
pile, are irreversible and a new "out-of-equilibrium" the- 
ory is required to describe them. In this paper we will 
describe the reversible states which are amenable to the 
ergodic hypothesis while the unconsolidated irreversible 
states will be addressed in a future work. 



III. VOLUME FUNCTION 

While it is always possible to measure the total volume 
of the system, it is unclear how to partition the space to 
associate a volume to each grain. Thus, the first step 
to study the V-ensemble is to find the volume function 
Wi associated to each particle that successfully tiles the 
system. This is analogous to the additive property of en- 
ergy in equilibrium statistical mechanics. This derivation 
is explained in Jamming I [43j and is the main theoretical 
result that leads to the solution of the partition function 
for granular matter. We refer to this paper for details. 
Here we just state the main two results, Eqs. ([S]) and 
([9]), and explain their significance for the solution of the 
statistical problem. 

Initial attempts at modelling W included, (a) a model 
volume function under mean-field approximation [3| , (b) 
the work of Ball and Blumenfeld [5|, |49|] in 2d and in 3d 
[50| . and (c) simpler versions in terms of the first co- 
ordination shell by Edwards [5l|. These definitions are 
problematic: (a) is not given in terms of the contact net- 
work, (c) is not additive, (b) and (c) are proportional to 
the coordination number of the balls contrary to expec- 
tation (see below). In Jamming I [4^ we have found an 
analytical form of the volume function in any dimensions 
and demonstrated that it is the Voronoi volume of the 
particle. 

The definition of a Voronoi cell is a convex polygon 
whose interior consists of all points closer to a given par- 
ticle than to any other (see Fig. [T]). Its formula in terms 
of particle positions for monodisperse spherical packings 
in 3d is 
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^ O^Ds, (8) 



where rij is the vector from the position of particle i to 
that of particle j, the average is over all the directions s 
forming an angle 9ij with as in Fig. [H and R is the 
radius of the grain. Wf defines the orientational Voronoi 
volume which is obtained without the integration over s. 
This formula has a simple interpretation depicted in Fig. 

m 




FIG. 1: The Voronoi volume is the light grey area (shown 
in 2d for simplicity). The limit of the Voronoi cell of particle 
i in the direction s is li{s) = rij /2 cos 6ij. Then the Voronoi 
volume is proportional to the integration of li(s)^ over s as in 
Eq. (El). 



The Voronoi construction is additive and successfully 
tiles the total volume. Prior to Eq. ([H]), there was no 
analytical formula to calculate the Voronoi volume in 
terms of the contact network rij . Equation ([8]) provides 
such a formula which allows theoretical analysis in the 
V-ensemble. However this microscopic version is difficult 
to incorporate into the partition function since it would 
necessitate a field theory. The next step is to develop a 
theory of volume fluctuations to coarse grain W™'' over a 
mesoscopic length scale and calculate an average volume 
function. The coarsening reduces the degrees of freedom 
to one variable, the coordination number of each grain, 
and defines a mesoscopic volume function which is more 
amenable to statistical calculations than Eq. ([5]). We 
find a (reduced) free volume function ^43i] : 



w{z) 



Va 



(9) 



valid for monodisperse hard spheres with volume Vg, 
where z is the geometrical coordination number (which 
is different from the mechanical coordination number, Z, 
see below). The average is over all the balls. 

The inverse relation with z in Eq. ([9]) is in gen- 
eral agreement with experiments (52l] . Many experi- 
ments have focused on the analysis of the system vol- 
ume and sin gle particle volume fluctuations in granular 
media [1, [ol. |52|. Is^ . Of particular importance to the 
present theory are the recent advances in X-ray tomog- 
raphy [s^l and confocal microscopy [23, [l^l which have 
revealed the detailed internal structure of jammed mat- 
ter allowing for the study of the free volume per particle. 
By partitioning the space with Voronoi diagrams, these 
studies show that Wi^"^ is distributed with wide tails 
[27I [51I [5^ . [53} . More importantly, the X-ray tomogra- 
phy experiments with 300,000 monodisperse hard spheres 
performed by the Aste group 52] find that the Voronoi 
volumes are inversely proportional, on average, to the co- 
ordination number of the particle, z, in reasonable agree- 
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ment with the prediction of the volume function Eq. ([9]). 
Such data is displayed in Fig. 6 in [52] where the local 
volume fraction defined as (j)^^ — 1 — Wi'""^ /Vg is plotted 
against the coordination number. From Eq. ^ we find: 



(10) 



1 z + 2a/3' 



which agrees relatively well with the shape of the curves 
displayed in Fig. 6 of [13] for different packing prepara- 
tions. It should be noted that for a more precise com- 
parison a coarse grained volume fraction and geometrical 
coordination number should be considered in Eq. pop . 
A detailed discussion is provided in Section IVlI Fl 

The free volume function decreases with z as expected 
since the more contacts per grain, the more jammed the 
particle is and the smaller the free volume associated 
with the grain. The coordination number z in Eq. ^ 
can be considered as a coarse-grained average associated 
with "quasiparticles" with mesoscopic free volume w. By 
analogy to the theory of quantum energy spectra, we re- 
gard each state described by Eq. Q as an assembly of 
"elementary excitations" which behave as independent 
quasiparticles. As a grain is jammed in the packing, it in- 
teracts with other grains. The role of this interaction be- 
tween grains is assumed in the calculation of the volume 
function ([9]) and it is implicit in the coarse-graining pro- 
cedure explained above. The quasiparticles can be con- 
sidered as particles in a self-consistent field of surround- 
ing jammed matter. In the presence of this field, the 
volume of the quasiparticles depends on the surround- 
ing particles as expressed in Eq. ([5]). The assembly of 
quasiparticles can be regarded as a set of non-interacting 
particles (when the number of elementary excitations is 
sufficiently low) and a single particle approximation can 
be used to solve the partition function as we will formu- 
late in Section IVTl 

While Eq. ([5]) is difficult to treat analytically, the ad- 
vantage of the mesoscopic Eq. ([9]) is that the partition 
function can be solved analytically since w depends on z 
only, instead of . The key result is the relation between 
the Voronoi volume and the coordination number which 
allow us to incorporate the volume function into a sta- 
tistical mechanics approach of jammed hard spheres, by 
using the constraint of mechanical stability as we show 
below. 



IV. DEFINITION OF JAMMING: ISOSTATIC 
CONJECTURE 

The definition of the constraint function 0jam is inti- 
mately related to the proper definition of a jammed state 
and should contain the minimum requirement of mechan- 
ical equilibrium. 

Distinguishing between metastable and mechanically 
stable packings that define the jammed state through the 
©jam function remains a problem under debate, related 



to the more fundamental question of whether or not a 
jammed packing is well-defined [s^. In practice, it is 
widely believed that the isostatic condition is necessary 
for a jammed disordered packing following the Alexander 
conjecture iSS, i57|, i58|] which was tested in several works 

It is well known that mechanical equilibrium imposes 
an average coordination number larger or equal than a 
minimum coordination where the number of force vari- 
ables equals the number of force and torque balance equa- 
tions [56, 53,(5^. The so-called isostatic condition. 

It is important to note that the derivation of the vol- 
ume function in Section IIIII implies nothing about the 
value of the contact forces; the volume function repre- 
sents the contribution arising purely from the geometry 
of the packing. Thus, the coordination number z appear- 
ing in Eq. ([9]) is the geometrical coordination number 
related to volume, which is different from the mechan- 
ical coordination number Z that counts the number of 
contacts per particle with non-zero force related to the 
isostatic condition and force network. 

In the case of frictionless spherical particles the iso- 
static condition is Z = 2d = 6 (in 3d), while the co- 
ordination in the case of infinitely rough particles (with 
interparticlc friction coefficient fi — > 00) is = d-l- 1 = 4, 
where d is the dimension of the system. Numerical sim- 
ulations and theoretical work suggest that at the jam- 
ming transition the system becomes exactly isostatic 
[13, 111 HH, m H, [13, Q . But no rigorous proof of this 
statement exists. In the following, we will use this iso- 
static conjecture to define the ensemble of jammed mat- 
ter. 

While isostaticity holds for perfectly smooth and in- 
finitely rough grains, the main problem is to extend it to 
finite frictional systems. For finite friction, the Coulomb 
condition takes the form of an inequality between the 
normal interparticlc contact force F„ and the tangential 
one Ft'. Ft < ^Fn and therefore no trivial solution to the 
minimum number of contacts can be obtained. 

The problem can be understood as an optimization 
of an outcome based on a set of constraints, i.e., mini- 
mizing a Hamiltonian of a system over a convex polyhe- 
dron specified by linear and non-negativity constraints. 
The isostatic condition can be augmented to indicate 
the number of extra equations for contacts satisfying the 
Coulomb condition, which is analogous to the number of 
redundant constraints in Maxwell constraint counting of 
rigidity percolation. This suggests that, for a finite value 
of /i, the original nonlinear problem can be mapped to 
a linear equation problem if we know how many extra 
equations should be added. This problem is treated with 
more details in a future paper. 

In the following, we present a rough approximation 
suggesting that the coordination number could be related 
to friction. We consider the following argument. Con- 
sider a set of spherical particles interacting via normal 
and tangential contact forces. These can be the stan- 
dard Hertz and Mindlin/ Coulomb forces of contact me- 
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TABLE I: Number of constraints and variables determining 
the isostatic condition for different systems of spfierical par- 
ticles. 



Friction 




Nt 


Ef 


Et 


^1 = 


^NZ 





dN 





fi finite 




i(d-l)7VZ/2(M) 


dN 


id(d-l)iV/i(M) 


/i = oo 


{NZ 


^{d-l)NZ 


dN 


id(d- l)iV 



chanics [22, 123, l2a| , respectively, see Section IVII Al We 
set N: number of particles, iV„: number of unknown 
normal forces, Nt: number of unknown tangential forces, 
Ef: number of force balance equations, Et: number of 
torque balance equations, Z — 2M/N: average coordina- 
tion number of the packing, where M is the total number 
of contacts, fi{fi): undetermined function of the friction 
coefficient /i such that 1 — is the fraction of spheres 
that can rotate freely (/i(0) = and /i(oo) = 1), and 
/2(/i): undetermined function of fi indicating the ratio of 
contacts satisfying Ft < /iF„, which satisfies /2(0) = 
and /2(oo) — 1. 

A packing is isostatic when 



Nn + Nt = Ef+ Et 



(11) 



The average coordination number at the isostatic point 
is then (see Table H]) : 



Z(^) = 2d 



l + l/2(rf-l)/i(A.) 
l + {d-l)f2i^i) ' 



(12) 



reducing to the known Z = 2d for frictionless particles 
and Z ^ d + 1 for infinitely rough particles. 

In what follows, we use the numerical fact that interpo- 
lating between the two isostatic limits, there exist pack- 
ings of finite fi with the coordination number smoothly 
varying between Z(/i = 0) = 6 and Z(/i —^ oo) — *■ 4 
[20| . Numerical simulations with packings described in 
Section IVlII corroborate this result and further show that 
the Z vs fi dependence is independent of the preparation 
protocol as obtained in our simulations (see Fig. [2]) . This 
result then generalizes the isostatic conditions from /i = 
and fi — oo to finite /i. 

It should be noted that we cannot rule out that other 
preparation protocols could give rise to other dependence 
of Z on jj,. However, we will see that the obtained phase 
diagram is given in terms of Z; the main prediction of 
the theory would still be valid irrespective of the par- 
ticular dependence Z{fi). That is, the theory does not 
assume anything about the relation between the inter- 
particle friction and Z. 

It is worth noting that other attempts to define the 
jammed state have been developed. A rigorous attempt 
is that of Torquato et al. who propose three categories 
of jamming [59|: locally jammed, collectively jammed 
and strictly jammed based on geometrical constraints. 




FIG. 2: Mechanical coordination number versus friction /i ob- 
tained in our numerical simulations explained in Section [VIII 
for different preparation protocols characterized by the initial 
volume fractions (jii indicated in the figure. The symbols and 
parameters used in these simulations are the same as in the 
plot of Fig. [11] 



A proper definition of jamming is intimately related to 
the existence of well-defined RCP and RLP limits [37l |. 
A well-defined RCP requires mathematical definitions of 
the jammed states and the concept of randomness [s^ . 
Indeed it was conjectured [5^ that the RCP concept 
could be replaced with the maximally random jammed 
(MRJ) conception in terms of an ensemble of order pa- 
rameters. 

However, the definition of [59] is based purely on geo- 
metrical considerations and therefore only valid for fric- 
tionless particles. Thus it is not suitable for granular 
materials with inter-particle frictional tangential forces; 
their configurational space is influenced by the mechanics 
of normal and shear forces. Other approaches to define 
a jammed state based on the potential energy landscape 
jl9l | fail for granular materials too since such a potential 
does not exist for frictional grains due to their inherent 
path-dependency. Thus, the definition of jammed state 
for granular materials must consider interparticle normal 
and tangential contact forces beyond geometry. In the 
companion paper j45| we elaborate on this problem. 



A. Geometrical and mechanical coordination 
number 

We have acknowledged a diff'erence between the ge- 
ometrical coordination number z in Eq. ([9]) and the 
mechanical coordination number Z which counts only 
the contacts with non-zero forces. Below we discuss the 
bounds of z. 
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Isostatic packing 



fi = Q z = 6 Z = 6 fi = oo 2 = 6 Z = 4 




FIG. 3: (a) Consider a frictionless packing at the isostatic 
limit with z — 6. In this case the isostatic condition implies 
also Z — 6 mechanical forces from the surrounding particles, 
(b) If we now switch on the tangential forces using the same 
packing as in (a) by setting fi oo, the particle requires only 
Z = 4 contacts to be rigid. Such a solution is guaranteed by 
the isostatic condition for /i — > oo. Thus, the particle still 
have z = 6 geometrical neighbors but only Z — 4 mechanical 
ones. 



ing to the volume function and the mechanical coordina- 
tion number, Z, measuring the contacts that carry forces 
only. This distinction is crucial to understand the sum 
over the states and the bounds in the partition function 
as explained below. 

These ideas are corroborated by numerical simulations 
in Section fVlIFl below. The packings along the vertical 
RCP line found in the simulations (see Fig. [TT|) have ap- 
proximately the same geometrical coordination number, 
z « 6. However, they differ in mechanical coordination 
number, going from the frictionless point Z — 6 to Z ^ 4 
as the friction coefficient is increased to infinity. 

We have established a lower bound of the geometri- 
cal coordination in Eq. (jl3p . The upper bound arises 
from considering the constraints in the positions of the 
rigid hard spheres. For hard spheres, the Nd positions 
of the particles are constrained by the Nz/2 geometrical 
constraints, \ fij\ — 2R, of rigidity. Thus, the number of 
contacts satisfies Nz/2 < Nd, and z is bounded by: 



Since some geometrical contacts may carry no force, 
then we have: 

Z <z. (13) 

To show this, imagine a packing of infinitely rough (/i 
oo) spheres with volume fraction close to 0.64. There 
must be z = 6 nearest neighbors around each particle 
on the average. However, the mechanical balance law 
requires only Z = 4 contacts per particle on average, 
implying that 2 contacts have zero force and do not con- 
tribute to the contact force network. 

Such a situation is possible as shown in Fig. [3] start- 
ing with the contact network of an isostatic packing of 
frictionless spheres having Z — % and all contacts car- 
rying forces (then z = 6 also as shown in Fig. [3^), we 
simply allow the existence of tangential forces between 
the particles and switch the friction coefficient to infinity. 
Subsequently, we solve the force and torque balance equa- 
tions again for this modified packing of infinitely rough 
spheres but same geometrical network, as shown in Fig. 

[Notice that the shear force is composed of an elastic 
Mindlin component plus the Coulomb condition deter- 
mined by /Lt, see Section IVH Al for details. Thus when 
/i — > oo, the elastic Mindlin component still remains]. 

The resulting packing is mechanically stable and is ob- 
tained by setting to zero the forces of two contacts per 
ball, on average, to satisfy the new force and torque bal- 
ance condition for the additional tangential force at the 
contact. Such a solution is guaranteed to exist due to the 
isostatic condition: at Z — 4 the number of equations 
equals the number of force variables. Despite mechanical 
equilibrium, giving Z = 4, there are still z = 6 geometri- 
cal contacts contributing to the volume function. 

Therefore, we identify two types of coordination num- 
ber: the geometrical coordination number, z, contribut- 



z < 2d. (14) 

Notice that this upper bound applies to the geometri- 
cal coordination, z and not to the mechanical one, Z, 
and it is valid for any system irrespective of the friction 
coefficient, from /i = — > cxd. 

Furthermore, in relation with the discussion of fric- 
tionless isostaticity, it is believed that above 2d, the sys- 
tem is partially crystallized. To increase the coordination 
number above 6, it is necessary to create partial crystal- 
lization in the packing, up to the point of full order of 
the FCC lattice with coordination number 12. Thus, 
by defining the upper bound at the frictionless isostatic 
limit we also exclude from the ensemble the partially 
crystalline packings. This is an important point, akin 
to mathematical tricks employed in replica approaches 
to glasses [28| . 

In conclusion, the mechanical coordination number, Z, 
ranges from 4 to 6 as a function of fi, and provides a lower 
bound to the geometrical coordination number, while the 
upper bound is 2d. A granular system is specified by the 
interparticle friction which determines the average me- 
chanical coordination at which the system is equilibrated, 
Z(/i). The possible microstates in the ensemble available 
for this system follow a Boltzmann distribution Eq. ([7]) 
for states satisfying the following bounds: 

Z{ii) <z<2d = 6. (15) 



V. DENSITY OF STATES 

According to the statistical mechanics of jammed mat- 
ter, the volume partition function Z is defined by Eq. ([7]). 
In the quasiparticle approximation we can write: 
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N (16) 

X 6jam 

i 

Considering N non-interacting quasiparticles with free 
volume w{z), the partition function can be written as: 

ZsX = ( / 9H ^'^"^ ©Jam dwf . (17) 

Here, g(vj) is the density of states for a given quasiparticle 
free volume. 

Since the mesoscopic w is directly related to z through 
Eq. ([9]), we change variables to the geometrical coordi- 
nation number in the partition function. The density of 
states for a single quasiparticle, g{w), is: 

g{w) = / P{w\z)g{z)dz, (18) 
Jz 

where P{w\z) is the conditional probability of a free vol- 
ume w for a given z, and g{z) is the density of states for 
a given z. Here, we have used the bounds in Eq. (fTS]! . 

The next step in the derivation is the calculation of 
g{z) which is developed in three steps: 

First, we consider that the packing of hard spheres is 
in a jammed configuration where there can be no collec- 
tive motion of any contacting subset of particles leading 
to unjamming when including the normal and tangential 
forces between the particles. This definition is an ex- 
tended version of the collectively jammed category pro- 
posed by Torquato et al. [59] that goes beyond the merely 
locally jammed configuration of packings, unstable to the 
motion of a single particle. While the degrees of freedom 
are continuous, the fact that the packing is collectively 
jammed implies that the jammed configurations in the 
volume space are not continuous. Otherwise there would 
be a continuous transformation in the position space that 
would unjam the system contradicting the fact that the 
packing is collectively jammed. 

Thus, we consider that the configuration space of 
jammed matter is discrete since we cannot change one 
configuration to another in a continuous way. Similar 
consideration of discreteness has been studies in [5^] . No- 
tice that the volume landscape could be continuous since 
we could change the volume as well, or, in the case of 
soft particle, we can deform them. Additionally, in the 
case of frictionless packings of soft particle, the energy of 
deformation is well-defined, and we can define the collec- 
tively jammed configuration as a minimum in the energy 
with definitive positive Hessian (or a zero order saddle) 
like in [l^. In this case, there could be no continuous 
transformation of the particle coordinates that brings one 
jammed state to the next, unless we deform the particles. 
Thus, the space is discrete in this case as well. 



Second, we call the dimension per particle of the con- 
figuration space as T) and consider that the distance 
between two jammed configurations is not broadly dis- 
tributed (meaning that the average distance is well- 
defined). We call the typical (average) distance between 
configurations in the volume space as hz, and therefore 
the number of configurations per particle is proportional 
to l/{hz)^ . The constant hz is very small [hz ^ 1) and 
plays the role of the Planck's constant, h, in quantum me- 
chanics. The situation is analogous to the discreteness of 
the configuration space imposed by the Heisenberg un- 
certainty principle in quantum mechanics. The formula 
is analogous to the factor h~'^ for the density of states in 
equilibrium statistical mechanics. While the degrees of 
freedom {pi, qi\ are continuous, the uncertainty principle 
imposes the discreteness (Ap, Aq) in the configurational 
space given by ApAq ~ h. Formally, the sum over the 
states Efj in equilibrium Ei, exp(— /Ji^iy) is replaced by, 

/ ex.]){-(3E{p,,qi))n^{dqdPt/h'^)- 

Third, we add z constraints per particle due to the fact 
that the particle is jammed by z contacts. Thus, there 
are Nz position constrains {\rij\ = 2R ) for a jammed 
state of hard spheres as compared to the unjammed "gas" 
state. Therefore, the number of degrees of freedom is 
reduced to P — z, and the number of configurations is 
then 1/(^12)^^^. Since the term \/{hz)'^ is a constant, 
it will not influence the average of the observables in the 
partition function (although it changes the value of the 
entropy, see Jamming HI fiSf). Therefore, we obtain 

g{z)^{hzY (19) 

with = \Ti{l/hz)- Physically, we expect /iz ^ 1, then 
g{z) is an exponentially rapid decreasing function with 
z. The most populated state is the highest volume at 
z = 4 while the least populated state is the ground state 
at z = 6. A constant could be added in (fTO|) but it has 
no consequence for the average of the observables in the 
partition function. However, it affect the value of the 
entropy. For the calculation of the entropy we consider 
that there is a single mesoscopic ground state and use 
g{z) = e-(^-2rf)Ac gg. 

We see that the negative of the geometrical coordi- 
nation number, — z, plays the role of a dimensionality 
or number of degrees of freedom for a packing, due to 
the extra position constrains of the contacting particles. 
The coordination z can be then considered as the num- 
ber of degrees of "frozen" per particle. Another way to 
understand Eq. is the following. In the case of a 

continuous phase space of configurations we would ob- 
tain ^'^g^^z)^ = 0. However, since the space of volume 
configurations is discrete as discussed above, the ratio 
^^lyi ~ hz- This imphes again Eq. P^ . 

The conditional probability P{w\z) depends on the w 
function, w — The average is taken over a certain 

mesoscopic length scale since the volume of a particle 
depends on the positions of the particles surrounding it. 
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Practically, such length scale is approximately of several 
particle diameters, w is a coarse-grained volume and 
independent of the microscopic partition of the particles, 
implying; 



P{w\z) ^5{w~2V3/z). 



(20) 



The meaning of Eq. (|^D|) is that we neglect the fluc- 
tuations in the coordination number due to the coarse 
graining procedure. A more general ensemble can be con- 
sidered where the fluctuations in the geometrical coordi- 
nation number are taken into account. We have extended 
our calculations to consider fluctuations in z and find a 
similar phase diagram as predicted by the present parti- 
tion function. This generalized z-ensemble will be treated 
in a future paper, where the boundaries of the phase di- 
agram are regarded as a second-order phase transition. 
The generalized z-ensemble allows for the calculation of 
the probability distribution of the coordination number, 
beyond the assumption of the delta function distribution 
in Eq. ((20)) . The resulting prediction of P{z) is in general 
agreement with simulations (see jS]). 

Substituting Eq. ^ and Eq. ^ into Eq. we 
find the isostatic partition function which is used in the 
remaining of this study: 



2^/3' 

zX 



dz. 



(21) 



VI. PHASE DIAGRAM 

Next, we obtain the equations of state to define the 
phase diagram of jammed matter by solving the parti- 
tion function. From Eq. (pij) . we calculate the ensemble 
average volume fraction cj) = (w + l)^^ — z/{z + 2'\/3) as: 



Ziso{X,Z)Jz Z + 2V3 



We start by investigating the limiting cases of zero and 
infinite compactivity. 

(a) In the limit of vanishing compactivity (X —^ 0), 
only the minimum volume or ground state at z = 6 con- 
tributes to the partition function. Then we obtain the 
ground state of jammed matter with a density: 



yRCP 



= cl>{X = 0,Z) = 



6 + 2\/3 



0.634, Z(^)G[4,6]. 

(23) 



The meaning of the subscript RCP in ([231) will become 
clear below. 

(b) In the limit of infinite compactivity {X oo), the 
Boltzmann factor exp[— 2V3/{zX)] 1, and the average 
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frictionless point 
(J-point) 
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FIG. 4; Phase diagram of jammed matter: Theory. 

Theoretical prediction of the statistical theory. All disordered 
packings lie within the yellow triangle demarcated by the RCP 
line, RLP line and G line. Lines of finite isocompactivity are 
in color. The grey area is the forbidden zone where no jammed 
packings can exist. 



in (I22p is taken over all the states with equal probability. 
We obtain: 



4>Rhp{Z) =(f>iX oo, Z) 



Ziso(oo, Z) Jz z + 2a/3 



: exp (z In hz) dz. 



(24) 



The constant determines the minimum volume in 
the phase space. We expect /i^ <C 1, such that the ex- 
ponential in Eq. ([M)) decays rapidly. Then the leading 
contribution to Eq. (I24p is from the highest volume at 
z — Z and therefore: 



z 



Z + 2\/3' 



Z(/i)e [4,6] 



(25) 



This dependence of the volume fraction on Z suggests 
using the (0, Z) plane to define the phase diagram of 
jammed matter as plotted in Fig. 31 The equations of 
state (1231) and (|25p are plotted in the (0, Z) plane in Fig. 
|4l providing two limits of the phase diagram. Since the 
mechanical coordination number is limited by 4 < Z < 6 
we have two more horizontal limits: The phase space 
is delimited from below by the minimum coordination 
Z = 4 for infinitely rough grains, denoted the granular- 
line or G-line in Fig. IH 

All mechanically stable disordered jammed packings 
lie within the confining limits of the phase diagram (in- 
dicated by the yellow zone in Fig. while the grey 
shaded area in Fig. 21 indicates the forbidden zone. For 
example, a packing of frictional hard spheres with Z = 5 
(corresponding to a granular material with interparticle 
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friction /i « 0.2 according to Fig. [J) cannot be equili- 
brated at volume fractions below (j) < 0rlp(-^ = 5) = 
5/(5 + 2V3) = 0.591 or above (f> > 0rcp = 0.634. It 
is worth noting that particular packings can exist in the 
forbidden zone; our contention is that they are zero mea- 
sure and therefore have zero probability of occurrence at 
the ensemble level. 

These results provide a statistical definition of the RLP 
and RCP limits: 

(i) Definition of RCP.— Stemming from the statis- 
tical mechanics approach, the RCP limit arises as the 
result of the relation (|23p , which gives the maximum vol- 
ume fraction of disordered packings under the mesoscopic 
framework. To the right of the RCP-line, packings exist 
only with some degree of order (for instance with crys- 
talline regions). The prediction, 
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^RCF — 



6 + 2V3' 



(26) 



or 63.4%, is valid for all friction coefficients and ap- 
proximates the experimental and numerical estimations 
psl [39I |40| . l4l| which find a close packing limit indepen- 
dent of friction in a narrow range around 0.64. 

Beyond the fact that 63-64% is commonly quoted as 
RCP for monodisperse hard spheres, we present a phys- 
ical interpretation of that value as the ground state of 
frictional hard spheres characterized by a given interpar- 
ticle friction coefficient. In this representation, as fi varies 
from to 00 and Z decreases from 6 to 4, the state of 
RCP changes accordingly while its volume fraction re- 
mains the same, given by Eq. (|26p . The present approach 
has led to an unexpected number of states that all lie in 
the RCP hue from Z = 6 to Z = 4 as depicted in Fig. gl 
suggesting that RCP is not a unique point in the phase 
diagram. 

This result may explain why experiments have found 
a consistent value of RCP for many distinct systems no 
matter the value of fi. The RCP can be found at the 
frictionless point with the highest possible coordination 
number Z — 6. Such a value of Z has been observed in 
the experiments of Brujic et al. on concentrated emul- 
sions near the jamming transition [60] and simulations of 
droplets. However, a value of 0rcp has not been deter- 
mined yet for monodisperse frictionless systems, in part 
due to the difficulties in synthesizing monodisperse emul- 
sions |27| . The more extensive evidence for a frictionless 
RCP appears from simulations, which has been exten- 
sively performed for hard and soft sphere systems. They 
find a common value of RCP for many different prepa- 
ration protocols [17|, ,1^, i3S] (there could exist some sub- 
tleties, see below). 

For frictional materials, experiments of Scott and Kil- 
gour |.39i] . and others display a nearly universal value of 
volume fraction for RCP consistent with the theoretical 
estimation. The experiments of Bernal [sj are interest- 
ing because they simultaneously measured the coordina- 
tion number and the volume fraction. Indeed, Mason, 
a postgraduate student of Bernal, took on the task of 



shaking glass balls in a sack and freezing the packing 
by pouring wax over the whole system. He would then 
carefully take the packing apart, ball by ball, painstak- 
ingly recordin g th e positions of contacts for each of over 
8000 particles [37j . Their result is a coordination number 
Z w 6 and a volume fraction <j) w 0.64 which corresponds 
to the prediction of the frictionless point in the phase di- 
agram. This may indicate that the particles used in the 
experiment of Bernal have a low friction coefficient. 

Another explanation, more plausible, is that the co- 
ordination number measured by Bernal is not the me- 
chanical one, but the geometrical one. Indeed, we may 
argue that the only coordination that can be measured in 
experiments of counting balls 'a la Bernal' is the geomet- 
rical one, since, one is never sure if the contacting balls 
were carrying a force or not. Bernal could only measure 
geometries, not forces. This is in addition to the uncer- 
tainty in the determination of the contacts with such a 
rudimentary method as pouring wax and then counting 
one by one the area of contact. The theory predicts that 
along the RCP-line, all the RCP states have geometrical 
z = 6, while Z ranges from 4 to 6. Thus, it is quite 
plausible that the coordination number in Bernal exper- 
iments is the geometrical one, z = 6, in agreement with 
the theory. 

It is worth noting that since the labor-intensive method 
patented half a century ago by Bernal, other groups have 
extracted data at the level of the constituent particles us- 
ing X-ray tomography 52]. Such experiments may give a 
clue to the relationship between coordination number and 
volume fraction. The experimental data, to date, seems 
in good agreement with the present theory. However, this 
method does not directly determine the contacts to verify 
whether the particles are touching or just very close to- 
gether. Furthermore, the method cannot measure forces 
so the distinction between geometrical and mechanical 
coordination is not possible to achieve in this experiment. 

An answer might be obtained using methods from bio- 
chemistry being developed at the moment. These meth- 
ods promise to provide the high resolution to determine 
the contact network with accuracy to develop an experi- 
mental understanding of the problem. 

An important prediction is that for frictionless sys- 
tems there is only one possible state at Z = 6. It is 
important to note that there is one state only at the 
mesoscopic level used in the theory. However, for a sin- 
gle mesoscopic state, we expect many microstates, which 
are averaged out in the mesoscopic theory of the vol- 
ume function. Thus, there could be more jammed states 
surrounding the frictionless point in the phase diagram. 
However, we expect that these states are clustered in a 
narrow region around the frictionless point. To access 
these microstates it would require different preparation 
protocols, analogous to the dependence of the glass tran- 
sition temperature on cooling rates in glasses [281 [soj . 

It is interesting to note that replica approaches to 
the jammed states [28^] predict many jammed isostatic 
states with different volume fractions for frictionless hard 
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spheres. The first question is whether these states in the 
force/energy ensemble are the same as the volume en- 
semble states that we treat in our approach. It should 
be noted that jamming in [28', 'so'l is obtained when the 
particles are rattling infinitely fast in their cages, in the 
limit of infinite number of collisions per unit time. That is 
when the dynamic pressure related to the momentum of 
the particles diverges. Second, from our point of view, an 
investigation of the fluctuations of the microstates as well 
as more general ensembles allowing for fluctuations in the 
coordination number could be considered. These studies, 
may reveal whether the frictionlcss point is unique or not. 
This point is investigated further in Section IVlI El 

(ii) Definition of RLP. — Equation of state ((25|) pro- 
vides the lowest volume fraction for a given Z and rep- 
resents a statistical definition of the RLP limit depicted 
by the RLP-linc in Fig. [H We predict that to the left 
of this line, packings are not mechanically stable or they 
are experimentally irreversible as discussed in [8, 9, 10]. 

A review of the literature indicates that there is no 
general consensus on the value of RLP as different es- 
timations have been reported ranging from 0.55 to 0.60 
[39L l4ll . [4^ , proposing that there is no clear definition of 
RLP limit. The phase diagram proposes a solution to this 
problem. Following the infinite compactivity RLP-line, 
the volume fraction of the RLP decreases with increasing 
friction from the frictionlcss point (0, Z) = (0.634, 6), to- 
wards the limit of infinitely rough hard spheres, Z — > 4. 
Indeed, experiments [s^ indicate that lower volume frac- 
tions are achieved for larger coefficient of friction. We 
predict the lowest volume fraction in the limit: fi 00, 
X — > 00 and Z ^ 4 (and — > 0) at 
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0.536. 



(27) 



Even though this is a theoretical limit, our results in- 
dicate that for /i > 1 this limit can be approximately 
achieved. 

The finding of a random loose packing bound is an 
interesting prediction of the present theory. The RLP 
limit has not been well investigated experimentally, and 
so far it was not certain whether this limit can or can- 
not be reached in real systems. The lowest stable volume 
fraction ever reported, 0.550 ±0.006, obtained by Onoda 
and and Liniger [4^ as the limit of vanishing gravity for 
spherical glass beads, is not far from the present predic- 
tion. 

The intersections of the RCP, RLP and the G-line iden- 
tify three interesting points in the {((>, Z{fi)) plane: 

(a) The frictionlcss point fi — 0, denoted J-point in 

at 

J= (0RCP,^(O)) = (0.634,6), 

corresponds to a system of compressed emulsions in the 
limit of small osmotic pressure as measured by J. Brujic 



(b) The lowest coordination number Z = i plotted as 
the G-line defines two associated points from the lowest 
volume fraction of loose packings at infinite compactivity, 
L-point, 

L=(0-["p,Z(oo)) = (0.536,4), 
to the zero compactivity state of close packing, C-point, 



C = 



ORCP 



Z(oo)) = (0.634,4). 



The full JCL triangle defines the isostatic plane where 
the frictional hard sphere packings reside. 

(in) Intermediate isocompactivity states. — For 
finite X, Eq. can be solved numerically. For each 

X, the function (l){X, Z) can be obtained and is plotted 
as each isocompactivity color line in Fig. |4l Between 
the two limits Eqs. (|23p and (PS)) . there are packings 
inside the yellow zone in Fig. 3] with finite compactivity, 
< X < 00. Since X controls the probability of each 
state, like in condensed matter through a Boltzmann-like 
factor in Eq. ([7]) , it characterizes the number of possible 
ways to rearrange a packing having a given volume and 
entropy, S. Thus, the limit of the most compact and 
least compact stable arrangements correspond to X ^ 
and X — > cxD, respectively. Between these limits, the 
compactivity determines the volume fraction from RCP 
to RLP. 

Dependence on and negative compactivity. 

Of interest is the dependence of our prediction on the 
"Planck constant" of jammed matter, /i^, that deter- 
mines the minimum volume in the phase space. The 
equation of state (|25p and the prediction of the mini- 
mum RLP, Eq. (j27p have been obtained by considering 
hz ^ but still nonzero. Thus, when X 00, the only 
state contributing to the volume partition function is the 
most populated at z = Z{p). 

The approximation /i^ <C 1 is a sensible one, since the 
discretization of the space is supposed to be very small. 
However, this constant remains a fitting parameter of the 
theory. Indeed in the simulations we will use a value of 
hz = e~^'^'^ in order to fit the theoretical values with the 
numerical ones for finite compactivity. This extremely 
small constant shifts the value of the minimum RLP a 
little bit to the right of the phase diagram from the pre- 
diction of Eq. (P7)) . Indeed, when we plot the phase 
diagram of Fig. |4] with a = 6"^*^° a slightly larger 
value of </'§Jl'p obtained as seen in Fig. [H In the un- 
physical limit of /iz ^ 1 we would obtain a minimum 
RLP value of ((^^^lp + '/'RCp)/2, although as said, this 
would correspond to an unphysical situation. 

However, it should be noted that this argument de- 
pends on the assumed exponential form of the density 
of states. Since at this point we do not have the exact 
form for the density of states, the results could change if 
other more accurate density is found to be valid. On the 
other hand the prediction of the ground state at RCP, 
Eq. remains unaffected by the density of state or 

the value oi h-,. 
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It is interesting to note that we can extend the com- 
pactivity to negative values and study the range X : 
0^ — !■ — oo. Indeed, for any value of h^, in the limit 
X ^ 0~ we obtain the lowest volume fraction of the 
prediction of Eq. ([25]) . Thus, properly speaking, the 
minimum value of RLP is defined in the limit of negative 
zero compactivity, and this limit is independent of the 
value of hz- The entropy has an interesting behavior in 
this regime which will be discussed later. 

It is important to realize that the region from X : 
0^ ^ oo disappears when —^ 0, thus the only mean- 
ingful limit is that of X ^ +oo (which equals both the 
limits X — oo and X 0^ when hz —^ 0). Thus, 
for any practical purpose, Eq. (I27p can be considered to 
be the lowest possible density predicted by the theory. 
However, notice that the shape of the iso-compactivity 
curves depends slightly on the value of hz. 

Since we always expect hz <^ I, it may not be neces- 
sary to use the negative compactivity states to describe 
RLP. Being sure, that it is still very difficult to accept the 
notion of compactivity in many quarters, even 20 years 
following its inception, and negative compaction may be 
an even more arbitrary concept at this point. However, 
we leave this interesting observation for further investi- 
gations where the entropy of jamming is treated in more 
details in [ist . 



A. Volume landscape of jammed matter 

Equation ([9]) plays the role of the Hamiltonian of the 
jammed system and the jammed configurations can be 
considered as the minima of ([9]). Inspired by the physics 
of glasses and supercooled liquids, we imagine a "volume 
landscape" , analogous to the energy landscape in glasses 
[6ll | , as a pictorial representation to describe the states of 
jammed matter. Each mesoscopic jammed state (deter- 
mined by the positions of the particles, denoted f;, and 
its corresponding volume w) is depicted as a point in Fig. 
[5l At the mesoscopic level, the volume landscape has dif- 
ferent levels of constant z, analogous to energy levels in 
Hamiltonian systems. 

The lowest volume corresponds to the FCC/HCP 
structure (with kissing number z = 12), as conjectured 
by Kepler [3^. Other lattice packings, such as the cubic 
lattice and tetrahedron lattice, have higher volume levels 
in this representation. Beyond these ordered states, the 
ensemble of disordered packings is identified within the 
yellow area in Fig. [5^, corresponding to a system with 
infinite friction. In this case the partition function is in- 
tegrated from z = 4 to 6, thus all the states are sampled 
in the configuration average as indicated by the arrow in 
Fig. [5^. When X ^ oo, all the states are sampled with 
equal probability, and, when X — > 0, the ground state 
is the most probable. As the compactivity is varied, the 
states along the G-line in the phase diagram result. For 
/i — > oo, the maximum volume, w{z = 4) — \/3/2 is at- 
tained for z = 4 when /i — > oo in analogy with the high 



energy states in a classical system. 

As we set the friction coefficient to a finite value, 
the available states in the volume are less, since the 
integration in the partition function is in the region 
Z{fJ.) < z < 6. This is indicated in Fig. [SJa for a generic 
Z{fi). Finally when the friction vanishes, we obtain only 
the ground state z = 6 as indicated in Fig. (St. 

For any value of /i, the lowest state is always at z = 6. 
This corresponds to the states exemplified by the RCP- 
line in the phase diagram, all of them with a geomet- 
rical coordination z = 6. Equation (j23p indicates that 
the RCP corresponds to the ground state of disordered 
jammed matter for a given friction which determines Z, 
while the RLP states are achieved for higher volume lev- 
els as indicated in Fig. [5l An important conclusion is the 
following: The states along the RCP line all have z = 6, 
independent of Z (and hz) which ranges from 4 to 6 as a 
function of the friction coefficient /i. The states along the 
RLP line all have z = Z (when hz 0). These predic- 
tions find good agreement with the numerically generated 
packings in Section IVlI Fl 

At fixed volume, the jammed states are separated by 
barriers of deformation energy, as depicted in Fig. 
These barriers can be understood as follows: so far 
we have treated the case of hard spheres considered as 
soft spheres (interacting via soft-potential such as Hertz- 
Mindlin forces) in the limit of vanishing deformation or 
infinite shear modulus. Indeed, the way to consider forces 
in granular materials is by considering the small tiny de- 
formation at the contact points and a given force law 
p^ . In a sense, deformable particles are needed when 
discussing realistic jammed states especially when con- 
sidering the problem of sound propagation and elastic 
behavior [H, [s^j ■ In the case of deformable particles the 
third axis in Fig. [SJi corresponds to the energy of defor- 
mation or the work done to go from one configuration to 
the next. This energy is not uniquely defined in terms 
of the particle coordinates; it depends on the path taken 
from one jammed state to the next. Thus, we empha- 
size that the energy in Fig. [5] is path-dependent. The 
only point where it becomes independent of the path is 
in the frictionless point. Besides this, the volume land- 
scape in the isostatic plane Figs. [5^-c is well defined and 
independent of the energy barriers and path dependent 
issues. 

It is important to note that the basins in Fig. [5] are 
not single states, but represent many microscopic states 
with different degrees of freedom r^, parameterized by 
a common value of z with a density of states g{z). The 
basins represent single states only at the mesoscopic level 
providing a mesoscopic view of the landscape of jammed 
states. This is an important distinction arising from the 
fact that the states defined by w{z) — 2\/3/z, Eq. ([9]), 
are coarse grained from the microscopic states defined 
by the microscopic Voronoi volume Eq. ([8]) in the meso- 
scopic calculations leading to (O as discussed in Jam- 
ming I psj . This fact has important implications for 
the present predictions which will be discussed in Section 
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FIG. 5: Schematic representation of the volume landscape of jammed matter (;fi,W). The multidimensional 
coordinate fi represents the degrees of freedom: the particle positions. Each dot represents a discrete mesoscopic jammed 
states at different z. It is important to note that for each meso state there are microscopic states with the same z. All 
disordered packings are in the yellow region (of an schematic shape) of the phase space which corresponds to the isostatic plane 
of hard spheres at the jamming transition where our calculations are performed. Other ordered packings have lower volume, 
such as the FCC. (a) We represent the case of /i = oo. The states represent those along the G-line in Fig. |4]as the compactivity 
varies from X = Q (ground state) to X ^ cxd at the RLP. The horizontal lines indicate packings at constant volume. The 
ground state of jammed matter for this friction coefficient has 2 = 6 and the highest volume states are found for z = 4. The 
arrow indicates the limits of integration in the partition function for this particular friction, (b) For another finite fj,, the space 
is delimited by above by a line of constant z — Z{fi). (c) For /x = only the ground state is available, giving rise to a single 
state, (d) The volume states states in (a,b,c) are separated by energy barriers represented by the third coordinate in the phase 
space. The energy barriers are path dependent due to friction between the particles. Nevertheless, the jammed configurations 
are well-defined in the isostatic plane, and the energy barriers represent the work done to go from one configuration to another. 
Only at the frictionless point, the energy barriers are path independent. 



IVII El The advantage of the volume landscape picture is 
that it allows visuahzation of the corresponding average 
over configurations that give rise to the macroscopic ob- 
servables of the jammed states. 



B. Equations of state 

Further statistical characterization of the jammed 
structures can be obtained through the calculation of 
the equations of state in the three-dimensional space 
{X, (/>, S), with S the entropy, as seen in Fig. [HI 

The entropy density, s = is obtained as: 

s{X,Z) = {w)/X + \iiZiUX,Z) (28) 

This equation is obtained in analogy with equilibrium 
statistical mechanics and it is analogous to the defini- 
tion of free energy: F = E ~ TS where F — — TlnZ is 



the free energy. We replace T ^ X , E ^ (w). There- 
fore, F ^ E - TS or S ^ {E - F)/T = E/T + hiZ 
is now s{X, Z) — (w) /X + lnZiso(^, Z), which is plot- 
ted as the equation of state in Fig. [51 For the calcu- 
lation of s we consider g{z) = /i^^"^'^ since for z = 2d 

there is only one mesoscopic state or g{z) ~ e with 
Zc = — 0.01 as used in Fig. [5] and [TlJ This cor- 

rection is irrelevant for average quantities like (/>(X, Z) 
but determines a constant factor in the entropy from the 
term In^iso. 

Each curve in the figure corresponds to a system with 
a different ^(/i). The projections S{X) and S'(0) in Fig. 
|6] characterize the nature of randomness in the packings. 
When comparing all the packings, the maximum entropy 
is at and X — > oo while the entropy is minimum for 

'/'RCP at X — > 0. Following the G-line in the phase dia- 
gram we obtain the entropy for infinitely rough spheres 
showing a larger entropy for the RLP than the RCP. 
The same conclusion is obtained for the other packings 
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FIG. 6: Predictions of the equation of state of jammed 
matter in the {X, (j), s) space. Each hne corresponds to a 
diflerent system with Z{jj,) as indicated. The projections in 
the (0, s) and {X, s) planes show that the RCP {X = 0) is less 
disordered than the RLP {X oo). The projection in the 
(X, (j)) plane resembles qualitatively the compaction curves of 
the experiments [EHiIlOl- 

at finite friction (4 < Z{ii) < 6). We conclude that the 
RLP states are more disordered than the RCP states. 
Approaching the frictionless J-point at Z = 6 the en- 
tropy vanishes. More precisely, it vanishes for a slightly 
smaller <p than (/)rcp of the order hz- Strictly speaking, 
the entropy diverges to — oo at 0rcp as S' ^ InX for 
any value of Z, in analogy with the classical equation of 
state, when we approach RCP to distances smaller than 
hz- However, this is an unphysical limit, as it would 
be considering distances in phase space smaller than the 
Planck constant. 

It is commonly believed that the RCP limit corre- 
sponds to a state with the highest number of configu- 
rations and therefore the highest entropy. However, here 
we show that the states with a higher compactivity have a 
higher entropy, corresponding to looser packings. Within 
a statistical mechanics framework of jammed matter, this 
result is a natural consequence and gives support to such 
an underlying statistical picture. A more detailed study 
of the entropy is performed in Jamming III [45| . 

The interpretation of the RCP as the ground state, 
X — > 0, with vanishing entropy (5* ^ and there- 
fore a unique state) warrants an elaboration. We no- 
tice that there exist packings above RCP all the way to 
'/'FCC = 0.74048, but these packings have some degree of 
order. These partially ordered packings do not appear 
in our theory because we treat only disordered packings 
characterized by the mesoscopic volume function which 
has been derived under the isotropic assumption. By do- 
ing so, we explicitly do not consider crystals or partially 



crystalline packings in the ensemble. This interprets the 
RCP in the context of the found third-law of thermo- 
dynamics. Our approach of neglecting the crystal state 
from the ensemble has analogies in replica treatment of 
glasses [62i] . 

Since at RCP the mesoscopic entropy vanishes, we 
identify this point as a granular "Kauzmann point" in 
analogy with the density (temperature) at which the con- 
figurational entropy of a colloidal (molecular) glass van- 
ishes at the ideal glass transition f61j. From this point 
the entropy increases monotonically with X, being max- 
imum for the RLP limit. We note that the microscopic 
states contribute still to the entropy of RCP giving rise 
to more states than predicted by the present mesoscopic 
approach. This case is discussed in more detail in 45]. 

The equation of state 4>{X) for different values of Z 
can be seen in the projection of Fig. [6l The volume 
fraction diminishes with increasing compactivity accord- 
ing to the theoretical picture of the phase diagram. The 
curves (f>{X) qualitatively resemble the reversible branch 
of compaction curves in the experiments of @ for shaken 
granular materials and oscillatory compression of grains 
[lOl suggesting a correspondence between X and shaking 
amplitude. The intention is that, different control param- 
eters in experiments could be related to a state variable, 
and therefore might help experimentalists to describe re- 
sults obtained under different protocols. 

For any value of Z , there is a common limit (f> 0rcp 
as X — !■ 0, indicating the constant volume fraction for 
all the RCP states. The singular nature of the friction- 
less J-point is revealed as the volume fraction remains 
constant for any value of X, explaining why this point 
is the confluence of the isocompactivity lines, including 
RCP and RLP. We conclude that at the frictionless J- 
point the compactivity does not play a role, at least at 
the mesoscopic level. 



C. Experimental realization of the phase diagram 

A reanalysis of the available experimental data tends 
to agree with the above theoretical predictions. How- 
ever, it would be desirable to perform more controlled 
experiments in light of the present results. As in other 
out of equilibrium systems, such as glasses, the inher- 
ent path-dependency of jammed matter materializes in 
the fact that different packing structures can be realized 
with different preparation protocols 0, H, H, 113 
volving tapping, fluidized beds, settling particles at dif- 
ferent speeds, acoustic perturbations or pressure waves 

Previous experiments [42l| and simulations ^251 find 
that lower volume fractions can be achieved for smaller 
settling speeds of the grains or slower compression (or 
quenching) rates during packing preparation. Colloidal 
glasses find a similar scenario [28l |3l|, with their glass 
transition temperature dependence on the quench rate of 
cooling, although due to different reasons. This raises 
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the question of whether the jamming point is unique or 
determined by the preparation protocol |31)] . 

The experiments performed by Onoda and Liniger with 
glass spheres in liquid of varying density to adjust condi- 
tions of buoyancy in the limit of vanishing gravity show 
that for larger values of gravity, the volume fraction de- 
creases. This indicates that settling speed of the par- 
ticles can determine the final volume fraction. Indeed, 
it was found numerically that the compression rate dur- 
ing preparation of the packings is a systematic way to 
obtain lower packing fractions, such that lower volume 
fractions can be achieved with quasi-static compression 
rates during preparation of the packings [2^ . 

These results provide a hint on how to generate the 
packings in the phase space and they are exploited in 
the following section treating the dynamical generation 
of the packings predicted in the phase diagram. Since Z 
is directly determined by /i, and the compactivity deter- 
mines what value of the volume fraction a packing has 
between the limits of the phase diagram, the main ques- 
tion is how to generate packings with different for a 
fixed /z to allow the exploration of the phase diagram. 



VII. NUMERICAL TESTS 

In this section we perform two different numerical 
tests: on the predictions of the theory and the assump- 
tions of the theory. The former are explained in Section 
IVII CI while the latter are elaborated in Section IVII Fl 
It is worth mentioning that while the predictions can be 
tested with packings prepared numerically and experi- 
mentally, the test of the assumptions of the theory is not 
so trivial. This is because the theory is based on the 
existence of quasiparticles which carry the information 
at a mesoscopic scale. Thus, in principle we cannot use 
real packings, such as computer generated packings or 
experimental ones, to measure the quasiparticles. The in- 
formation obtained from those packings already contains 
ensemble averages through the Boltzmann distribution 
and density of states. That is, it is in principle not possi- 
ble to isolate the behavior of quasiparticles from the real 
measures rendering difficult to properly test some of the 
assumptions of the theory at the more basic mesoscopic 
scale. A proper way to test the quasiparticles is by gener- 
ating fully random packings, a study developed in Jam- 
ming [i^. Nevertheless, in Section IVlI Fl we attempt to 
perform such a test, especially since we can easily obtain 
the packings at the RLP line or the infinite compactiv- 
ity limit. These packings are fully random obtained as 
flat averages in the ensemble without the corresponding 
Boltzmann factor and may contain direct information on 
the mesoscopic fluctuations. 

The existence of the theoretically inferred jammed 
states opens such predictions to experimental and com- 
putational investigation. We numerically test the predic- 
tions of the phase diagram by preparing monodisperse 
packings of Hertz-Mindlin [23, l2g| spheres with friction 



coefficient \i at the jammin g tr ansition using methods 
previously developed [13, [Ulil [111 . 

We test the theoretical predictions and show how to 
dynamically generate all the packings in the phase space 
of conflgurations through different preparation protocols. 
Although diverse states are predicted by the theory, they 
may not be easily accessible by experimentation due to 
their low probability of occurrence. For instance, we will 
flnd that the packings close to the C-point (high volume 
fraction, high friction, low coordination) are the most 
difflcult to obtain. 

The advantage of our theoretical framework is that 
it systematically classifles the different packings into a 
coherent picture of the phase diagram. In the follow- 
ing we use different preparation protocols to generate all 
the phase space of jamming. In particular we provide a 
scheme to reproduce the RCP and RLP-lines amenable 
to experimental tests. 

We achieve different packing states by compressing a 
system from an initial volume fraction with a com- 
pression rate F in a medium of viscosity (damping) r\ 
where the particles are dispersed. The system is defined 
by the friction coefficient /i which sets according to 
Fig. [21 While the simulations are not realistic (no grav- 
ity, boundaries, or realistic protocol is employed), they 
provide a way to test the main predictions of the theory. 
The final state (0, Z) is achieved by the system for every 
((/>i, F, ry, /i) at the jamming transition of vanishing stress 
with a method explained next. 

It should be noted that other experimental protocols 
could be also adapted to the exploration of the phase 
diagram, including (a) Gentle tapping with servo mecha- 
nisms that adjust the system at a specified pressure [gij, 
(b) Gentle tapping with external oscillatory perturba- 
tions, (c) Settling of grains under gravity in a variety of 
liquids with viscosity ?/, (d) fiuidized beds. 

Relation vifith hard sphere simulations. — The 
present algorithm finds analogies with recent attempts 
to describe jamming using ideas coming from the the- 
ory of mean-field spin glasses and optimization problems 
[28l [sij . This is an interesting situation since it has been 
shown that in the case of hard spheres the system al- 
ways crystallizes unless an infinite quench is applied [ssj . 
However, our soft sphere simulations do not produce ap- 
preciable crystallization in the packing. 

This situation can be understood as following: In order 
to avoid the (partial) crystallization in practice, a fast 
compression is needed. A simple dimensionality analysis 
tells that the quenching rate F oc G^^^, where G is the 
shear modulus of the grains (see Fig. [7]). This analysis 
allows for comparison of our soft ball simulations with 
the hard sphere simulations used in other studies [ssj . 

One could imagine that there exists a typical quench 
rate, depending on the stiffness of the grains, as a division 
of the ordered and disordered (jammed) phases. It is 
interesting to note that at the hard-sphere limit, G — > oo, 
all the packings are in the ordered phase except when 
F 00, where (/> « 0.64. This explains the behavior of 
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FIG. 7: Dimensional analysis indicates that for finite shear 
modulus G of the particles, the system can crystallize or ran- 
domize according to the quenched rate F. This argument al- 
lows for a comparison with hard-sphere simulations that are 
done in the limit G — » oo. In this case only F — > oo avoids 
partial crystallization. 

the hard sphere packings found in previous works [55| . 
For the soft ball case, the situation is more subtle since 
a certain range of compression rate is acceptable for the 
acquirement of disordered packings. It is in this range 
that we perform simulations. 

According to this picture, in the theoretical limit of 
infinitely slow quench rates, the system should in princi- 
ple crystallize for any G. In granular matter, this limit 
is almost impossible to achieve, since it is always found 
that the system is random and far from the crystal struc- 
ture. This is in part because, in the absence of tem- 
perature bath, the granular system is always stuck in a 
jammed configuration far from the crystal one. The con- 
tention is that, for a well defined set of quenches, the sys- 
tem does not crystallize and the study of the disordered 
phase ensues. Indeed, numerical and experimental evi- 
dence points to the validity of this assertion. The present 
theory neglects the existence of the ordered phase and 
assumes an ensemble of the disordered states. A more 
general theory will account for the ordered states as well. 
But this is a very difficult problem in jamming: a full 
theory that accounts for RCP and crystallization. Some 
ideas are presented in Jamming IV [46l |. From the sim- 
ulation point of view, we can say that our soft sphere 
approach does not produce crystallization while the hard 
sphere approach does not produce randomness for any 
finite compression rate. 

A. Molecular dynamics simulations of grains 

We prepare static packings of spherical grains interact- 
ing via elastic forces and Coulomb friction (see [ISjIl^ for 
more details). The system size ranges from TV = 1,024 
to = 10, 000 particles. Two spherical grains in contact 



at positions xi and X2 and with radius R interact with a 
Hertz normal repulsive force |25|] 

= ^ kr.Ry^6^'\ (29) 
and an incremental Mindlin tangential force [2^ 

AFf = A:t(i?(5)i/2As^ (30) 

Here the normal overlap is (5 = (1/2) [2i?— |xi — a?2|] > 0. 
The normal force acts only in compression, Fn = Q when 
(5 < 0. The variable s is defined such that the relative 
shear displacement between the two grain centers is 2s. 
The prefactors fc„ = 4G/(1 - v) and h = 8G/(2 - u) 
are defined in terms of the shear modulus G and the 
Poisson's ratio v of the material from which the grains 
are made. We use G = 29 GPa and v ~ 0.2 typical 
values for spherical glass beads and we use i? = 5 x 10~^ 
m and the density of the particles, p = 2 x 10^ kg/m^. 
Viscous dissipative forces are added at the global level 
affecting the total velocity of each particle through a term 
—72? in the equation of motion, where 7 is the damping 
coefficient related to the viscosity of the medium rj = 
7/(67r_R). Sliding friction is also considered: 

Ft < fiF„. (31) 

That is, when Ft exceeds the Coulomb threshold, iiFn, 
the grains slide and Ft = fiFn, where fj, is the static 
friction coefficient between the spheres. We measure the 
time in units of Iq = Ry^ p/G, the compression rate in 
units of Fq = 5.9t(^^ and the viscosity in units of 770 = 
8.2i?^p/to. The dynamics follows integration of Newton's 
equation for translations and rotations. 

B. Preparation protocol: Packings at the jamming 
transition with the split algorithm 

The critical volume fraction at the jamming transition, 
4>c (which is called in the theoretical part of this pa- 
per, we will use both notations for convenience) will be 
identified by the "split" algorithm as explained in [2^ . 
allowing one to obtain packings at the critical density 
of jamming with arbitrary precision. Preliminary results 
indicate that lower values of (j) are obtained for slow com- 
pression [l^l . Thus we expect that by varying F or 77 the 
phase space will be explored. Also, the initial plays 
an important role in achieving packings close to the RCP 
line. 

The preparation protocol consists of first preparing a 
gas of non-interacting particles at an initial volume frac- 
tion (f>i in a periodically repeated cubic box. The particles 
do not interact and therefore the stress in the system is 
CT = and Z = 0. 

To generate a random configuration with friction at 
volume fraction ~ 0.64 is very difficult. Therefore to 
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FIG. 8: Schematic representation of the protocol used to dy- 
namically achieved the packings in the phase diagram, (a) 
The region near the C-point (high friction, high volume frac- 
tion) is the most difficult to access. Here, we prepare an un- 
jammed sample with no friction at (=0.62 in the diagram) 
very close to RCP, and then we switch friction to compress to 
the required frictional jammed state. The packing follow the 
paths indicated in the figure for two friction coefficients. On 
the other hand, the RLP-line is easily obtained by compress- 
ing unjammed packings with low (j>i (=0.52 in the diagram) 
at (b) a given compression rate Fi for a given friction coef- 
ficient, (c) There is a small dependence on F, the faster the 
compression the deeper we enter in the phase diagram. 



achieve any volume fraction in this initial stage, we ini- 
tially work with a frictionless system. We first generate 
a very dilute (unjammed gas of non interacting particles) 
grain configuration randomly, usually with a volume frac- 
tion <j)f) « 0.30 ~ 0.36. Then we apply an extremely slow 
isotropic compression without friction on this dilute con- 
figuration until the system reaches another unjammed 
configuration at higher density, (see Fig. |5^). The 
reason to add an extremely slow isotropic compression to 
the dilute configuration is to avoid involving any kinetic 
energy to keep the system as random as possible. After 
obtaining this unjammed state with initial volume frac- 
tion (f>i, we restore friction and a compression is applied 
with a compression rate V until a given volume fraction 
(j)i. Then the compression is stopped and the system is 
allowed to relax to mechanical equilibrium by following 
Newton's equations without further compression. 

After the compression to , two things can occur (see 
Fig. [9|): 

(a) The system jams: If the system is above the jam- 
ming point (f>i > (j>c, then the stress will decrease and 
ultimately stabilize to a finite nonzero value, meaning 
that the pressure of the system remains unchanged (usu- 
ally Act < IQ-^ Pa, red fine, FigO over a large period of 
time (usually ^ 10^ MD steps). The coordination num- 
ber usually has a first initial decrease, but if the system 
is jammed it will also stabilize at a constant value above 
the isostatic minimal number (inset Fig. [5]). 



FIG. 9: Time evolution of stress (the pressure in the system) 
for two packings simulated as explained in the text. The 
solid red line represents a packing with <j>i > 0c and dotted 
black line represents a packing with <j>2 < 0c, where 0i = 
(j>2 + 2 X 10~*. The inset shows the time evolution of the 
coordination number. 

(b) The system is not jammed: here the stress and 
the coordination number will relax to zero. This fact 
is illustrated in Fig. [9l If the packing has 0i > 4'c, 
it stabilizes at a non-zero pressure above the jamming 
transition, but the pressure decreases very quickly to zero 
(the system is not jammed) if 02 < 4'c, even though 0i 
and 02 differ only by 2 x 10"*. 

We identify the exact volume fraction of the jamming 
transition 0c, as follows [l!]: A search procedure consist- 
ing of several cycles is applied such that in each cycle we 
fix the lower and upper boundaries of 0c, see Fig. [TOl 
The difference between the boundaries gets smaller as 
the cycles proceed, meaning that 0c is fixed with higher 
and higher precision. We start from the packing of high 
volume fraction 0i > 0c and generate a series of pack- 
ings with step-decreasing volume fractions until the first 
packing with zero pressure is observed, which has a vol- 
ume fraction = 0i — A0. Thus, 0c is bounded between 
01 — A0 and 0i. Then we test = 0i — A0/2. If 

= 01 — A0/2 is stable, 0c is between 0i — A0 and 

01 — A0/2. If = 01 — A0/2 is unstable, 0c is between 
01 — A0/2 and 0i. Therefore, in this cycle, we reduce 
the region where 0c possibly lies in from A0 to A0/2. 
If we carry out this cycle for n times, we improve the 
precision to A0„ = A0/2". In our simulations cycling 
ceases when A0„ gets below 2 x 10"'* and n = 12. A 
similar algorithm was employed in (2^ to study the ap- 
proach to the jamming transition by preparing packings 
at a finite pressure. In the present work we are inter- 
ested in jammed packings at vanishing pressure, right at 
the jamming transition in the isostatic plane defined for 
all friction coefficients. 

It is important to determine whether the packings are 
jammed in the sense that they are not only mechani- 
cally stable but also they are stable under perturbations. 
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FIG. 10: Flow chart for searching procedure of critical volume 
fraction (j>c. 



Our numerical protocols assure that the system is at 
least locally jammed since each particle is in mechani- 
cal equilibrium [s^ . To test if the system is collectively 
jammed is more involved. For frictionless systems, where 
tangential forces are removed, we use the Hertz energy 
C^hoitz = ^ knR^/'^S^/'^ to test whether the Hessian of 
the jammed configurations is positive [l^. We find that 
the frictionless configurations have positive Hessian indi- 
cating that they are collectively jammed. However, this 
method is not useful when considering frictional systems. 
The energy of deformation depends on the path taken to 
deform the system and cannot be defined uniquely. In 
this case, a numerical test of the stability of the packings 
applies a small random velocity to each jammed particle. 
We find that the packings are stable to small pertur- 
bations consisting of external forces of the order of 0.1 
times the value of the average force, indicating that the 
packings may be collectively jammed. 

To test for more strict jammed conditions involves 
studying the stability under boundary deformations. We 
have tested that our packings are stable under the most 
common strain deformations to isotropic packings by per- 
forming a uniaxial compression test, a simple shear and 
a pure shear test. 

A simple shear test implies a strain deformation 
Aei2 — Ae2i 7^ 0, while the rest of the strain compo- 
nents tij remain constant. A pure shear test is done 
with Aeii = — Ae22 7^ and a uniaxial compression test 
along the 1-direction is performed by keeping the strain 
constant in Ae22 = Aeaa = 0, and Aen ^ 0. Here, the 
strain e^-, is determined from the imposed dimensions of 
the unit cell. For example, en — AL/Lq where AL is 
the infinitesimal change in the 11 direction and Lq is the 
size of the reference state. 

In all cases the packings are stable under strain per- 
turbations indicative of the existence of a finite shear 
modulus (from the shear and compression tests) and a 
finite bulk modulus (from the compression test) of the 



jammed packings. A full investigation of the elasticity 
of the jamming phase diagram is left for future studies. 
It suffices to state that the numerically found states can 
be considered to be mechanically stable jammed states 
amenable to a statistical mechanical description in the 
sense of the reversible branch of compaction in the ex- 
periments of [1, [IS [13. 



C. Results: phase diagram 

We repeat the above procedure 10 times with differ- 
ent random initial configurations to get a better aver- 
age of 0c- Figure [TT] shows the results for a system of 
N — 1024. All of the packings shown in the present 
work as well as the code to generate them are available at 
http://jamlab.org. Each data point corresponds to a sin- 
gle set of {4>i^ r, 77, /i) and is averaged over these 10 reahza- 
tions. We consider 0.40 < <j)i < 0.63, 10"^ < F < 10"^, 
10^^ < T] < 10~^, and < < 00. The error of (j)c 
obtained over the 10 realizations as shown in Fig. [11] is 
usually 5 x 10"''', larger than the precision of the split 
algorithm (2 x 10^"'). In general, all numerically gen- 
erated jammed states lie approximately within the pre- 
dicted bounds of the phase diagram. 

The plot in Fig. [TTk explores the dependence of the 
packings {4>,Z) on the initial state 4>i. (In the plot, 
refers to the (j)c obtained in a split algorithm). In Fig. 
[TT]we plot our results for a fixed quench rate F — 10"'' 
and damping coefficient 77 = 10~^ (except for the last 
orange curve on the right with 77 = 10~^) and for differ- 
ent initial states ranging from left to right (see Fig. [Tl] 
for details) Each color here corresponds to a particular 
4>i (0.40 at the left to 0.63 at the right), while each data 
point along a curve corresponds to a system prepared at 
a different friction from /i = 0atZ«6to/i— s-cxaat 
Z Ri 4. For instance, the states joined by the dashed lines 
in Fig. [TTk are for systems with the same /i. While these 
lines are approximately horizontal, indicating an approx- 
imate independence of Z on the preparation protocol (pi, 
we observe deviations specially around the RLP line for 
intermediate values of fj,. Despite this, the Z{fj,) depen- 
dence is approximately valid as seen in Fig. [2l which is 
plotted using the results of Fig. [TTk . 

We find that the packings prepared from the larger ini- 
tial densities (f>i closely reproduce the RCP line of zero 
compactivity at (/)r,cp, confirming the prediction that 
RCP extend along the vertical line from the J-point to 
the C-point. 

We find that the packings along the RCP line have 
equal geometrical coordination number 2: « 6 but differ in 
their mechanical one from Z = 6 to « 4, in agreement 
with theory (see Section [VII Fp . These states are then 
identified with the ground state of jammed matter for a 
given /i as depicted in the volume landscape picture of 
Fig. [3 

In the limit of small initial densities (see curves for 
— 0.40, 0.53 and 0.55 in Fig. [TTk ) we reproduce ap- 
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FIG. 11: Phase diagram of jammed matter: Simula- 
tions. Numerical simulations demonstrate how to dynam- 
ically access the theoretically found states. The numerical 
protocol is parameterized by {(j>i,T,ri, fi). (a) The plot shows 
the dependence of the final jammed states {(j>, Z) on 0i for a 
fix r = 10~^ and r; = 10"'^ (except for the orange ► curve 
which is for rj = 10"*) from left to right 4>i — 0.40 (black ■), 
0.53 (red •), 0.55 (violet A), 0.57 (blue T), 0.59 (green ♦), 0.61 
(pink m) and 0.63 (orange ►). Each equal-color line set repre- 
sents a different (j>i and the dashed lines join systems with the 
same friction fi. Solid lines represent the theoretical results 
with hz = 6"^°" or Zc = 0.01 for different compactivities mea- 
sured in units of W~^Vg. From left to right X = oo (black 
solid line), 1.62 (blue solid line), 1.38 (green solid line), 1.16 
(pink solid line) and 0.88 (orange solid line). The error bars 
correspond to the s.d. over 10 realizations of the packings, 
(b) The plot focuses on the dependence of (0, Z) on (F, rj) for 
two different (jii. Solid black symbols are for (pi = 0.40 and 
(r,T;) = (lO-^lO-^) (black ■) and (10-^ 10"*) (black .). 
Open red symbols are for (jji = 0.63 and (F, 77) — (10~^, 10""^) 
(red □) and (10"^ 10"") (red o). 



proximately the predictions of the RLP-line. These pack- 
ings follow the theoretical prediction for infinite com- 
pactivity except for deviations (less than 5%) in the co- 
ordination number for packings close to the lower value 
of Z — A, where the numerical curve seems to flatten out 
deviating from the theory. Numerically, we find the low- 
est volume fraction at 4>'^£p = 0.539 ± 0.003, close to the 
theoretical prediction, (fi^^ = 0.536. We notice that the 
theoretical lines in Fig. [Tl]are obtained for hz = 6"^"°, 
thus the minimum RLP deviates slightly from the the- 
oretical prediction in the limit /i^ — )■ as discussed in 
Section ED 

Packings with the lowest Z correspond to infinitely 
rough spheres. The deviation of coordination number 
between theory and simulation (specially at low volume 
fraction) could be from the system not achieving an iso- 
static state at infinite friction. We were not able to gen- 
erate packings with exactly Z = A but our algorithm 
generates packings very close to this isostatic packings 
at Z w 4.2 for /i = 10**. In general, we find that while 
there are many states along the RLP line, the barriers 
of these states decrease as the coordination number de- 
creases towards Z ~ 4, i.e. when the friction increases. 
Thus, the states at the lower left part of the phase dia- 
gram are very difficult to equilibrate. 

Besides the states delimiting the phase space, we gen- 
erate other packings with intermediate values of 4>i = 
0.57, 0.59, 0.61 as shown in Fig. [TTk . Interestingly, we 
find that these states (all obtained for fixed F = 10~^ 
and 77 = 10"'^) closely follow the predicted lines of iso- 
compactivity as indicated in the figure. We find that the 
simulations from <j>i = 0.57, 0.59, 0.61 correspond to com- 
pactivities X — 1.62,1.38, and 1.16, respectively (mea- 
sured in units of 10^"^ Vg). 

The constant weakly affects the finite compactiv- 
ity lines. We find that = 6"^°" provides the best fit 
to the data for finite isocompactivity lines in Fig. [TTk . 
Thus, with reasonable approximation and for this partic- 
ular protocol, we identify the density of the initial state, 
4>i with the compactivity of the packing, providing a way 
to prepare a packing with a desired compactivity. We 
note, however, that other values of hz produce approxi- 
mately the same phase diagram boundaries (i.e., the RLP 
and RCP lines, as long as <C 1) but do not fit the iso- 
compactivity lines within as accurately. Thus, the iden- 
tification of the compactivity with (f)i remains dependent 
on this particular value of hz used to fit the theory. 

We also test the dependence on the state of the pack- 
ings with the compression rate F and viscosity rj. Both 
parameters should have similar effects since they slow 
down the dynamics of the grains. Figures [SId and[8j; rep- 
resent typical preparations. In general, we reproduce the 
RLP line for slow quenches or for large viscosities as seen 
in Fig. [ITb (for — 0.40 in black). In this case, the 
grains are allowed to slide and develop large transverse 
displacements and Mindlin forces with the concomitant 
low (j) and large compactivity. Therefore we find a pre- 
dominance of the Mindlin forces over the normal Herztian 
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forces as a characteristic of the lower volume fractions 
of RLP, having a pronounced path dependent structure 
with large compactivity in the packing. When the com- 
pression rate is increased or the damping is reduced we 
find packings with higher volume fractions as indicated 
in Fig. [Tib . Indeed, analogous protocol to generate these 
packings would be to change the damping coefficient by 
immersing the particles in liquids of different viscosities. 

It is also interesting to investigate the packings gener- 
ated by fixed friction and varying quenched rate, since 
this situation corresponds to a given system. For a given 
fi we find approximately a common value of Z{^) inde- 
pendent of r and (f)i, as discussed in Fig. [21 indicating 
that the curve Z{^) is roughly independent of the proto- 
col. Thus, a fixed corresponds to a horizontal line in 
the phase diagram as indicated in Fig[TTk by the dashed 
lines. 

Figure [8| shows a summary of the numerical protocols 
we use to generate the packings in the phase diagram. 
The "easiest" packings to generate are along the RLP line 
(Fig. [HJs, c) while the most difficult are deep in the phase 
diagram with low coordination number (large friction) 
and high volume fraction, that is near the C-point (Fig. 
[HK). In general, we always start with a dilute unjammed 
(non-interacting gas) sample without friction and com- 
pressed slowly until another dilute sample not jammed 
at a volume fraction 0j. If (pi is around or smaller than 
0.55 then the RLP line is obtained. If we decrease by 
3 orders of magnitudes the compression rate using this 
initial density, we move inside the phase diagram, though 
not substantially, as indicated schematically in the figure. 

Therefore, in order to access the inner states around 
the C-point at high friction, we took another route. 
We prepared a frictionless dilute unjammed sample at a 
higher </),i, which can be as close as possible to ~ 0.63 since 
frictionless spheres only jam at RCP. Then we switched 
on friction and got the packings inside the phase dia- 
gram as indicated in Fig. [8^. We realize that switching 
on and off friction is impossible in experiments. However, 
a possible way to experimentally realize the path to the 
C-point is by resetting the path dependency to zero along 
the preparation protocol: this path could be reproduced 
by stabilizing a frictional packing with fast compressions 
and resetting the interparticle forces by gently tapping. 
This will simulated the lost of path dependence that we 
simulate by the initial compression without friction and 
may allow to reach the packings around the C-point. 



D. Compactimeter for granular matter and the 
ABC experiment 

At this point we do not have a theoretical explanation 
for why packings with the same initial unjammed state (pi 
end up having the same compactivity when jammed after 
compression. It is important to note that this finding is 
based on the particular value we use for h^, such that a 
change in gives rise to different lines of isocompactivity 




FIG. 12: Measuring the temperature of sand. Pictorial rep- 
resentation. 

which may not fit the simulation results of Fig. [TTk so 
closely. Nevertheless, the close fit between theory and 
simulation deserves an explanation. 

We may conjecture that <f>i determines a type of dis- 
order quenched in the initial configuration that leads to 
systems with the same compactivity but different volume 
fractions and coordination numbers, evidenced by our re- 
sults. We can use this empirical result to control the com- 
pactivity of the packing, at least for this particular pro- 
tocol, defining a "granular compactimeter" . Laboratory 
measurements of compactivity usually involve indirect 
measures through Fluctuation-Dissipation relations be- 
tween fiuctuations and response functions [l,[l,[i3| • How- 
ever, there is no simple thermometer (or "compactime- 
ter" ) to measure this observable directly in a packing. 
Furthermore, there is no simple way to prepare a packing 
with a desire compactivity. Our results can be used to, 
at a minimum, define packings with equal compactivity. 
The empirical identification of X with (pi promotes the 
possibility of controlling X within this particular proto- 
col, an important step in any thermodynamical analysis. 

The advance of a granular thermodynamics crucially 
depends on the invention of a thermometer that can eas- 
ily measure the compactivity, as in Fig. 1121 A zero-th law 
of thermodynamics presents serious challenges in granu- 
lar materials since there is no straight forward way to 
define a compactivity bath or reservoir. 

The present theory predicts the dependence of the vol- 
ume fraction, (p{X), on the compactivity, allowing for the 
development of a compactimeter: Such a device consists 
of a known granular medium, having a given /i and, there- 
fore, a given coordination number. The compactimeter 
is inserted in a granular system and both, system and 
compactimeter, are gently shaken to allow for equilibra- 
tion at a common compactivity (the compactimeter has 
a flexible membrane that allows volume transfer between 
both systems in contact). The volume of the grain in 
the compactimeter will adjust accordingly and the com- 
pactivity can be read directly from the scale attached to 
it through the equation of state, Ar(0), provided by the 
theory. 
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FIG. 13; Sketch depicting a possible ABC experiment to test 
the zero-th law of granular thermodynamics and the validity 
of the phase diagram for finite compactivities. For interme- 
diate frictions, the A and B systems should follow the arrows 
to C when put into contact and allow to interchange volume 
by gently shaking. 



The foundation of the above compactimeter relies on 
the validity of a zero-th law for granular thermodynamics 
that determines the equilibration at a given compactiv- 
ity of a system: Two systems in contact should mechan- 
ically equilibrate at the same compactivity. The test of 
the zero-th law is difficult to facilitate without theoretical 
guidance. Thus, it will be useful to perform the follow- 
ing test, either numerically or experimentally to test not 
only the idea of equilibration but also the possibility to 
describe granular matter under the V-ensemble. Such a 
test has been labeled the ABC experiment as proposed 
by Edwards (sij . which is now possible to perform using 
the phase diagram. 

Two packings are numerically prepared at points A and 
B in the Fig. [13] for different /i, /ia < Mb- The systems 
are equilibrated at different compactivities < Xb 
by using different They are then put into contact 
through a flexible membrane (or by simply putting the 
particles at the surface in contact) and allowed to me- 
chanically equilibrate by gently shaking the system. Nu- 
merically, this can be simulated with the same compres- 
sion methods as explained in Section IVIII The total 
volume is kept fixed and the volumes of the subsystems 
should change accordingly, V — Va + Vb- If the AB sys- 
tem equilibrates at the same compactivity, then it will 
follow the trajectory depicted in Fig. [T3] towards the 
state C of equilibration along the isocompactivity line. 
There is no need to measure the compactivity of the final 
packing. By just measuring the final volumes Va and Vb 
(or their volume fractions) will suffice. This will provide 
an important test for the V-ensemble and compactivity 
to describe jammed matter, as well as the validity of the 
phase diagram. 

Another intriguing possibility would be to mix a fric- 



tionless A packing at the J-point with an infinite friction 
B packing at the L-point, which can be achieved numer- 
ically. This AB system would provide the maximum dif- 
ference between the isolated packings, allowing study of 
the zero-th law with more accuracy. According to theory, 
the frictionless packing is independent of X . Therefore, 
it should remain at the J-point. The fi oo packing 
should, in principle, equilibrate at any X or volume frac- 
tion along the G-line. Most probably, it would stay put at 
the L-point since the A and B packings are joined by the 
isocompactivity X oo line already. In any case, the re- 
sulting AB packing should fall inside the phase diagram 
if the theory is correct. If it does not equilibrate (falling 
outside the phase diagram, we note that the two packings 
are at the borderline), it could indicate the breakdown 
of the approach. Two conclusions could be reached from 
the failure of such a test. Either the V-ensemble is wrong 
or the approximations of the V-ensemble theory arc not 
correct. In the latter, more sophisticated theories should 
be developed. If the former case occurs, then a different 
ensemble will have to be considered beyond isostaticity. 



E. Microscopic states versus mesoscopic states 

The notion that RCP arises for X = 0, implying that 
there arc no fluctuations at RCP, deserves some discus- 
sion. Furthermore, it is important to discuss the implica- 
tions of the related coincidence of RCP and RLP at the 
frictionless point. That is, the fact that for frictionless 
particles there is a unique state of jamming and therefore 
the compactivity is irrelevant when /i = 0. 

First of all, we need to emphasize that these results 
should be interpreted at the mesoscopic level. Thus, 
there is one mesoscopic state at RCP and only the meso- 
scopic entropy vanishes at this point. These results pro- 
vides an unambiguous interpretation of RCP only at the 
meso level. As discussed above, for a single mesoscopic 
state, there are many microscopic states average out in 
the calculation of the average free volume function w{z) 
which is used in the partition function Eq. pT|) (see also 
Jamming I [43]). Therefore, we expect that these mi- 
croscopic states will contribute to the entropy of RCP 
(see Jamming III for a discussion) as well as giving rise 
to other states perhaps with different volume fraction. 
Below we elaborate on the existence of microscopic and 
mesoscopic packing states. This point is better discussed 
in terms of the volume landscape of Fig. [S) all the 
jammed states are degenerated around the mesoscopic 
ground state with the coordination number z — 6. These 
states have slightly different volume fractions, leading to 
microscopic fluctuations which are coarse-grained in the 
mesoscopic theory. 

Next, we investigate more carefully the microscopic 
states numerically by focusing on the dependence of the 
volume fraction at the frictionless J-point J = (0.634, 6) 
on the system size and on different preparation protocols 
by changing the initial volume fraction of the unjammed 
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FIG. 14: Investigation of the microscopic states of jamming 
for a frictionless system. We prepare packings with different 
number of particles ranging from = 256 up to = 10, 000 
particles. From A'^ — 256 to = 4096 we perform 10 realiza- 
tions of the packings to calculate the critical volume fraction 
at jamming. The error bars are s.d. over these realizations. 
For A'^ — 10* we have only one realization (for large A'^, the 
system is self-averaging) . We mainly focus on the dependence 
of (j> on the density of the initial unjammed state, <j>i. In gen- 
eral when (j>i increases close to (RCP) 0.63, the final density 
increases as shown. The result is valid independent on the 
size of the packings. There is a slight dependence on the 
compression rate F or the damping -q as shown. 



state (f)i, the compression rate F and rj. 

Figure [14] shows the dependence of the volume frac- 
tion at jamming for /i = as a function of the system 
size from N = 256 to N =10,000. Judging from the cur- 
rent numerical results, a series of packings with a small 
range of volume fraction from « 0.632 to « 0.642, but 
roughly the same coordination number 6, could be gen- 
erated by changing the (jji small dependence on the 
compression rate F is also observed in the third digit in 
(j) but more systematic studies are needed to corroborate 
this result for larger N). The lowest jamming density, 
(p « 0.632 is obtained for the smaller initial density of 
the unjammed initial packing (j>i = 0.40 while the largest 
density at jamming (j) ^ 0.642 is obtained when the initial 
density is close to RCP, (pi — 0.63, as shown. The max- 
imal volume fraction (p ~ 0.642 for frictionless packings 
is suggested as an upper bound for all the volume frac- 
tions under different preparations. We wish to point out 
that more systematic size dependence studies are been 
developed to obtain more precise values. 

It is worth noting that our numerical simulations indi- 
cate that these packings exist also in the frictional region 
of the phase diagram up to /i — > oo. We conclude that as 
the initial volume fraction of the unjammed frictionless 
state, <pi, gets close to RCP, the final volume fraction cpc 
for frictionless case starts to deviate from the theoretical 
prediction. It seems there is no typical value of volume 



fraction at J-point, which fits all different preparations, 
although the deviations are small. 

This result can be interpreted as a manifestation of 
the microscopic states neglected at the mesoscopic level. 
This raises an interesting analogy with recent work of 
Zamponi and Parisi [2^ where the jammed states are 
considered as infinite pressure glassy states obtained after 
a fast compression from a liquid state at finite tempera- 
ture. They propose that a range of different volume frac- 
tions of jamming can be achieved according to the initial 
state from where the quench is started, with the maxi- 
mum possible jammed density obtained when the initial 
density is the Kauzmann point of the ideal glass. In their 
representation (see Fig. 



4 in larXiv:0808. 2180^1 [28|) a 
range of initial densities in the liquid phase (p G [ipd, 'Pk], 
where ipd is the density where many metastable states 
first appear in the liquid phase as suggested by mean 
field model picture of the glass transition, and fx is the 
Kauzmann density of the ideal glass transition, gives rise 
to a final density of jamming, (fj{ip), obtained after a fast 
quench from ip. Two limiting densities of jamming are 
predicted: the minimum at ipj{ipd) — fth and the max- 
imum volume fraction obtained after a quench from the 
Kauzmann point, ipj{ipK) = fo- 

The situation is somehow analogous to our simula- 
tions and it is tempting then to identify tpth ~ 0.632 
and ipo ~ 0.642 with our numerical values. However, we 
consider this comparison premature. First, we need to 
investigate more deeply the size dependence of our re- 
sults to test whether the two lines in Fig. [14] could meet 
at oo. We note that they slightly overlap within 

error bars over 10 realizations. Second, it is still a ques- 
tion whether the states at infinite pressure obtained us- 
ing fast quenches from a liquid state "a la Zamponi" are 
equivalent to the jammed state obtained in our approach 
which are mechanically equilibrated states. Jamming in 
the hard-sphere approach of jl^ means particles moving 
infinitely fast in their cages; a very different concept to 
jamming in our soft particle approach, where particles 
are mechanically equilibrated. 

Nevertheless, an analogous critique could be said to our 
approach using soft balls since jamming is achieved when 
pressure or the average force vanishes. Besides these dif- 
ferences, our volume ensemble approach might still dif- 
fer from the energy/force ensemble used in the work of 
Zamponi-Parisi. Finally, we notice that the state with 
final density w 0.642 is difficult to obtain since we have 
to start from a very high density of the unjammed state, 
thus rendering the lower bounds as the most probable 
states to obtain in experiments, explaining their experi- 
mental ubiquity. Overall, these problems raise interesting 
questions to be investigated in the future. 

How can the microscopic states be modeled beyond the 
mesoscopic theory described in the present work? We are 
currently investigating this question. It seems that an ex- 
act description of the probability distribution of the mi- 
croscopic volumes is possible (at least up to a prescribed 
coordination shell of particles) providing a framework to 
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investigate the microscopic states. Furthermore, our the- 
ory suggests that the elementary units or building blocks 
of jammed matter are the quasiparticles of fix coordina- 
tion number. An exhaustive investigation of these local 
packings may provide clues to the above questions. 

The distinction between microscopic and mesoscopic 
implies an intrinsic difference between the current meso- 
scopic theory, Eq. (PT|) . and the full Edwards' statistics 
given by Eq. ([7]). A very tentative scenario is as follows. 
We could separate the entire entropy S to two parts as 
S — 5*111030 + •S'lnicro respectively, in terms of the differ- 
ent length scales. For the frictionless packings, we have 
seen that S'meso vanishes and then S = S'micro and the 
microscopic compactivity Xmicro = — is equivalent 
to the Edwards compactivity: X^dw = ■ Thus, one 

compactivity X^icm = A^Edw is needed to describe the 
jamming point in the frictionless case. 

The frictional case is more complicated since both 
mesoscopic and microscopic entropies play important 
roles. The fact that there exist packings beyond the pre- 
dicted mesoscopic volume fraction of RCP with S'meso = 
implies that there must be another compactivity to de- 
scribe at least these packings. 

Thus, we suggest a two-compactivity scenario based 
on the separation of microscopic and mesoscopic 
length-scale, which is different from the Edwards/one- 
compactivity scenario valid for frictionless grains. The 
FDT then provides a formalism to test the fluctuations 
and measure both compactivities which is written accord- 
ing to different scales: 



(^^nicro) ^micro 



incso/ mcso 



d{V) 
d{V) 



(32) 



and the total volume fluctuation 



xi 



OAii 



-){V). (33) 



Experimentally, the compactivities could be accessible 
by the intensity of tapping. We expect that the mi- 
croscopic states have smaller "volume" barriers than the 
mesoscopic states, thus they should be accessible via very 
gentle tapping that do not change the global Z, while 
stronger tapping would give the mesoscopic states. 



F. Geometrical versus mechanical coordination 
number: measuring the behavior of quasiparticles 



In Section IVII CI we have tested the predictions of the 
theory. In this Section we attempt to test some of the 
assumptions of the theory, specially those related to the 
quasiparticles of w{z) and the bounds in z and its distri- 
bution. A quasiparticle is a mathematical entity that is, 
in principle, impossible to isolate from the structure of 



a real packing prepared either with MD computer sim- 
ulations or experiments, rendering the measurement of 
their properties from real packings difficult. This is be- 
cause real packings already contain the ensemble average, 
and we can not dissociate the ensemble, with the Boltz- 
mann distribution of states, compactivity and density of 
states, from the isolated statistics of a quasiparticle. 

Quasiparticles can be tested by preparing random 
packings with prescribed geometrical coordination num- 
ber. However, below, we argue that we may obtain 
approximate properties from real packings in the limit 
X oo and ^ based on the following observa- 
tions: First, when X ^ oo the Boltzmann factor is one, 
and the average in the partition function is flat over the 
conflgurations up to the density of states. Furthermore, 
the density of states is found to be a very fast decay- 
ing function of z. These two results implies that along 
the RLP-line of real packings (either numerical or ex- 
perimental) we can obtain fully random packings which 
could reveal the properties of the quasiparticles. Under 
this approximation we can then test some of the approx- 
imations of the theory. However, the conclusions from 
this section remain approximate. 

We have proposed that the bounds of the geometrical 
coordination number are Z < z < 6 and below we test 
these bounds. We find numerically that the quasiparti- 
cles size is of the order of two particle diameters. Once z 
is averaged over this range, we find that z ranges approxi- 
mately between (Z, 6). In fact, we find even more narrow 
distribution than theoretically expected. This is due to 
the fact that we define z without the microscopic fluctua- 
tions but including the mesoscopic fluctuations. Once we 
average over its neighbors, we not only remove the mi- 
croscopic fluctuations but also, partially, the mesoscopic 
ones. In the limiting case of no mesoscopic fluctuations, 
when X = along the RCP line, we flnd a very narrow 
distribution at z w 6 for any value of Z and ^, after 
averaging z over a mesoscopic region of two particle di- 
ameters. The distribution is even narrower when z is 
coarse-grained over a region of four particle diameters. 

The X-ray tomography experiments of Aste (Fig. 
6a) reveal the trend predicted by Eq. ^ between the in- 
verse of the Voronoi volume and the number of neighbors 
of a set of Voronoi cells with similar volumes. However, it 
is evident from the flgure that such number of neighbors 
is spanning a range of values between 4 and ~ 10. The 
reason why the Aste's group 52] found a wider range of 
z is possibly due to the fact they did not consider a range 
of coarse-graining. 

A question arise as how to identify the geometrical 
from the mechanical coordination number, even numeri- 
cally. Strictly speaking, this distinction is not possible to 
materialize in real packings. Again, we need to generate 
random packings from a geometrical point of view, with- 
out resorting to forces, and then study their distribution 
in real packings where forces are taken into account. Nev- 
ertheless, we follow the following approximate scheme to 
try to obtain information of the quasiparticles from real 
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FIG. 15: gz{/S.r) of packings with various friction coefficient 
/J, along RCP line. We set 2R^1. 
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FIG. 16: gz{Ar) of packings with various friction coefficient 
fj, along RLP line. 

packings. 

Even in numerical simulations, there is always a round 
off error associated to measurements. Thus, two particles 
that may not be in contact (giving rise to a zero force) 
may be close enough to be considered as contributing to 
the geometrical coordination. Indeed, it is known that 
g{r) has a singularity [65], g{r) ^ (r — .5)~'''^, implying 
that there are many particles almost touching. Following 
this consideration, we introduce a modified radial distri- 
bution function (RDF) gz{r) in order to approximately 
identify z and Z from real packings: 

(34) 

where R is the radius of particle, N is the number of 
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FIG. 17: PDF of the coarse-grained coordination number {z}i 
for packings with various friction coefficient fi along the RCP 
line, (a) Ar = and « = 2; (b) Ar = 0.04 and I = 2; (c) 
Ar = 0.04 and I = 4. 



particles, r^j is the distance of two particle's centers, 
rij = \ri — rj\, and 6 is the Heaviside step function. 
The RDF describes the average value of the number of 
grains in contact with a virtual particle of radius r, and 
the factor of /r^ is the ratio of a real sphere's area and 
the virtual one's. 

gz (r) measures the number of balls with their volume 
intersecting the surface of a sphere of radius r measured 
from the center of a given ball. When r = i? in ([34]) 
we obtain the mechanical coordination number while the 
geometrical one is obtained for a small value Ar = ^ 
for which we distinctly find a signature from computer 
simulations, unambiguously defining it at Ar = 0.04. 

Figures [TSl and [TBI plot the gz{Ar) of packings with var- 
ious friction coefficient /i on the isostatic plane along the 
RCP and RLP lines respectively. Following the definition 
of Eq. ([M]) . gz, with Ar = 0, should be directly equal 
to the mechanical coordination number, Z, and should 
range from 4 to 6 along both RCP and RLP (if hz -^1) 
lines which is confirmed by our numerical simulations in 
Figs. [Island [TBI respectively. 

Furthermore, as shown in the figures, we find that 
gz{Ar) along the RCP line, increases slightly as Ar in- 
creases, and finally reach the same value of gz at Ar — 
0.04 as shown in Fig. [151 We identify the geometrical 
coordination number as z = 5^(0.04) under the accuracy 
of the simulations. Therefore, we conclude that all RCP 
states have approximately the same geometrical coordi- 
nation number, z « 6, in agreement with the theory. In 
terms of the volume landscape of Fig. [SI the states along 
the RCP line are ground states with different friction at 
X = 0. Thus they should all have z = 6. We notice 
that gz{Ar = 0.04) = 6.65 > 6, which may result, firstly, 
from the approximate nature of our measurements of z 
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FIG. 18: PDF of the coarse-grained coordination number (z); 
for packings with various friction coefficient fi along the RLP 
line, (a) Ar = and Z = 2; (b) Ar = 0.04 and I = 2. 
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FIG. 19: Summary of the theoretical findings regarding the 
range of z and Z along the different iso-z, iso-Z, and iso-X 
lines and the J, C, and L-points in the phase diagram. 



from real packings since we are not measuring the quasi- 
particles directly, and, secondly, from the increasing of 
the coordination number slightly away from the critical 
point as the local volume fraction cj) increases slightly as 
Z — Zc ^ ((/> — (/>c)''- We note that Z is also slightly 
increased from the isostatic point. 

Along the RLP line. Fig. [161 we find that the geomet- 
rical coordination number as extracted from f;z(0.04) is 
very close to the mechanical one. Since the RLP line is 
at X ^ oo and /iz <C 1, all the states along RLP have 
z K, Z a& we move along the line varying the friction co- 
efficient. Thus, the numerical results confirm the theory, 
and further confirm that the value of is very small. 



Next, we study the coarse-grained coordination num- 
ber (z)i as defined as. 



Ni N, 



R 



:-i e 



R 



--1 



> R 



(35) 

where Ni is the number of the particles inside a coarse- 
grained spherical range with a radius of /. Figures 1171 
and [THl plot the PDF of {z)i for all the packings along 
the RCP and RLP lines, respectively. The distributions 
show a narrow shape, and the average value gives the 
value of Qzir), i.e., 



{z)i ^ 9z{r). 



(36) 



We find that all the P{{z)i) along the RCP line coincide 
at Ar = 0.04 as shown in Fig. [T7b and [T7b . demon- 
strating that all the states have approximately the same 
average value of z (as discussed above), as well as the 
same distribution. The distribution gets narrower as the 
coarse-graining parameter / increase from Z = 2 to Z = 4 
as shown in Fig. [T7b . On the other hand Z changes as 
the friction changes following the isostatic condition, as 
seen in Fig. 117b .. This in in good agreement with the 
theory. 

Along the RLP line, the situation is analogous. As 
discussed above, z ^ Z . This is clearly demonstrated in 
Fig. [TSb which should be compared with the analogous 
Fig. [TTb for the RCP line. Figure [T8k also shows how 
Z changes with friction in the same way as in Fig. 117b . 
More importantly, we find that the bounds of (z); are well 
approximated between the bounds proposed by the the- 
ory {Z,6). The slight shift towards higher coordinations 
observed in Fig. 118b is due to the approximate nature 
of our measurements from real packings and the small 
increment in volume fraction from the critical point of 
jamming, as discussed above. Overall, we conclude that 
along the RCP line, all states have z « 6, with Z chang- 
ing from 6 to 4 as friction is increased. On the other 
hand, along the RLP line, we find that z « Z. These 
two results confirm the theory very well. Figure [1^] sum- 
marizes all these finding indicating the ranges of z and Z 
along the iso-X lines of RCP and RLP, the iso-z vertical 
lines, and the iso-Z horizontal lines in the phase diagram. 



VIII. OUTLOOK 

We have presented several results in the disordered 
sphere packing problem and here will discuss their sig- 
nificance. We start with the results that we believe are 
more robust and may survive higher scrutiny, and follow 
with the results that require further study plus a discus- 
sion on how to improve the theory. 

The final result of the phase diagram, the characteri- 
zation of RCP and RLP using the compactivity, the role 
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of friction in the determination of the hmiting packings 
are interesting results, ahiiost independent of the theory. 
The distinction between the geometrical and mechanical 
coordination number is also important, allowing for the 
implementation of averages in the partition function. On 
the mathematical side, the microscopic volume function 
and the formula obtained for the Voronoi cell in Jam- 
ming I (43} are important and exact results which are 
the foundation of all preceding work presented herein. 

While the predictions Eqs. IJ) and ^ of RCP and 
RLP approximate very well the known values, they are 
obviously not exact expressions, since they are based on 
several approximations of the statistical theory of vol- 
umes. Thus, the implementation of the mesoscopic vol- 
ume function and the probability distribution of volumes 
on which these results are based are, perhaps, the points 
where improvements may be necessary. While the 
formula reproduces the average volume fraction surpris- 
ingly well, it is based on approximations to calculate the 
probability distribution of volumes, P>(c), as explained 
in Jamming I [i^. More work is needed to obtain bet- 
ter approximations to such a distribution to capture not 
only the mean value of the Voronoi volume but higher 
moments as well. A full discussion is performed in Jam- 
ming IV (46j . The present theory shows the way to im- 
prove upon Jamming I and obtain an exact formulation of 
P> (c) , at least to a prescribed value of the coordination 
shell. Such a formulation is presently being developed 
and is based on a systematic study of quasiparticles with 
fixed geometrical coordination number; it produces good 
fits of P{W-) for the case of 2d, while 3d calculations are 
under way. 

Other approximations that require additional atten- 
tion are the considerations leading to the density of states 
g{z). While the main results of RCP and RLP are nearly 
independent of the particular form of 5(2), other quan- 
tities such as the entropy and the equations of state for 
finite nonzero compactivity are more sensitive to it. By 
a direct measurement of the entropy, g{z) can be calcu- 
lated and compared to the simple exponential form used 
in the present study. 

Furthermore, we have explicitly removed crystal states 
from our calculations by considering the upper bound of 
z = 6. A more general theory is needed then to investi- 
gate the transition from the RCP to the FCC as crystals 
starts to form in the system. Such a theory will include 
not only disordered states but also partially crystallized 
states and would give X = at FCC following the full 
Edwards partition function, Eq. ([7]). Whether this dis- 
order/order transition is critical or not is an interesting 
question which is under investigation. The present the- 
ory provides some ideas on how to proceed. The study of 
the distributions of local packings with a fix geometrical 
coordination number may provide some ideas to under- 
stand this problem. 

The distribution of coordination number, assumed to 
be a delta function, is another approximation to be im- 
proved upon. Once again, this may not affect the pre- 



dicted RCP and RLP but deeper studies require also a 
distribution of z. A more general ensemble is being con- 
sidered that predicts the distribution of z in good agree- 
ment with the existing data (see Jamming IV [46*1). 

Finally, the extension to study the microstates requires 
more investigation. The identification of the elementary 
units as quasiparticles of fixed coordination number is 
interesting and could eventually lead to more precise so- 
lutions of RCP and RLP. 



IX. CONCLUSIONS 

In conclusion, using Edwards statistical mechanics we 
have elucidated some aspects of RLP and RCP in the 
disordered spherical packing problem. The numerical re- 
sults suggest a way to experimentally test the existence 
of the predicted packings. By allowing the grains to set- 
tle in liquids of varying density, the speed of the particles 
can be varied and a systematic exploration of the jam- 
ming phase diagram can be achieved. 

Beyond the elucidation of some questions in the sphere 
packing problem, other problems can now be addressed 
systematically from the point of view of what we have 
learned from the phase diagram. This includes the inves- 
tigation of the criticality of the jamming transition from 
frictionless to frictional systems by extending the phase 
space to (0, Z, a), and generating packings away from the 
isostatic plane by isotropically compressing the packings 
generated in the plane of Fig. 31 Previous studies have 
focused on the frictionless point and on fixed friction. 
According to our results it is important to sample all 
the phase space for different friction and compactivities 
according to the preparation protocols explained above. 

Previous results find power law scaling 
where the stress vanishes as cr ~ ((/) — (jic)"- A possible 
scenario is depicted in the extended phase diagram of 
Fig. [201 in the 3-dimensional space (0, Z, a), which could 
be tested numerically and experimentally. We note in 
passing that, all the results in the present study refer to 
the isostatic plane at the jamming transition cr = or 
hard sphere limit. We expect that, as we compress pack- 
ings of soft spheres, different planes of cr = constant 7^ 
will be obtained as shown in Fig. [201 These compressed 
states are described not only by the compactivity X but 
also the angoricity A through the full partition function 
in the VF ensemble, Eq. ([4]). Our preliminary results 
indicate that there seems to be evidence of criticality as 
the isostatic plane is approached when cr — > 0, although 
with exponents dependent on friction as well as the value 
of the compactivity. 

Additional topics to be addressed by means of the pre- 
sented approach include characterization of jamming in 
the phase space of configurations, the problem of elastic- 
ity and Green's function and the study of the pair cor- 
relation function g{r), the volume distribution P(Wf), 
Voronoi volume, P{W^™), and coordination number 
P{z) , which are within the reach of the present approach. 
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FIG. 20: Sketch depicting a possible extension of the phase 
diagram of Fig. Uin the {(j), Z) isostatic plane at the jamming 
transition to the full phase diagram in the (0, Z, a) space to 
investigate the criticality of the jamming transition for all 
packings with any friction coefficient. 

Other systems under consideration are two-dimensional 
systems as well as the mean field case of infinite dimen- 
sion. These are interesting cases, since there are many 
questions debated as to what is the more dense state 
in higher dimension. Extensions to other systems such 
as anisotropic particles are also within the reach of the 



present approach. 

Finally, a complete characterization of the VF- 
ensemble can now be performed through the measure- 
ment of other quantities as a function of the volume 
fraction, such as the force distribution P(F„,Ft), and 
equations of state through the angoricity in the ((/>, Z, a) 
space. Some of these quantities can be obtained analyti- 
cally from the ensembles allowing for interesting predic- 
tions. 

The phase diagram introduced here serves as a begin- 
ning to understand how random packings fill space in 3d. 
The comparative advantage of the present approach over 
extensive work done in the past, is in the classification of 
all packings through X, Z and </> in the theoretical phase 
diagram from where these studies could be systemati- 
cally performed. This classification guides the search for 
indications of jamming from a systematic point of view, 
through the exploration of all jammed states from /i = 
to /i — > oo. 
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